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This dissertation discusses solutions to the nonlinear Klein-Gordon (nlKG) 
equation with symmetric and asymmetric double- well potentials, focusing on the 
collapse and collision of bubbles and critical phenomena found therein. 

A new method is presented that allows the solution of massive field equa- 
tions on a (relatively) small static grid. A coordinate transformation is used that 
transforms typical flatspace coordinates to coordinates that move outward (near the 
outer boundary) at nearly the speed of light. The outgoing radiation is compressed 
to nearly the Nyquist limit of the grid where it is quenched by dissipation. The 
method is implemented successfully in both spherically symmetric and axisymmet- 
ric codes. 



VI 



The new method is first used in a code to explore spherically symmetric 
bubble collapse. New resonant oscillon solutions are found within the solution space 
of the nlKG model with a symmetric double- well potential (SDWP). A time-scaling 
relation is found to exist for the lifetime of each resonance. The resonant solutions 
are also obtained independently using a set of ordinary differential equations derived 
from a non-radiative periodic ansatz. The method is also applied to the nlKG model 
with an asymmetric double- well potential (ADWP); the threshold of expanding 
bubble formation is investigated and a time-scaling law is shown to exist. 

The method is then used in an axisymmetric code to simulate bubble col- 
lisions. A technique for boosting arbitrary spherically symmetric finite difference 
solutions is presented and used to generate initial data for the collisions. The 2D pa- 
rameter space of bubble width versus collision velocity is explored and the threshold 
of expanding bubble formation is again considered. On the threshold, there exists 
a time-scaling law with critical exponent similar to the spherically symmetric case. 

Lastly, resonant oscillon solutions are constructed using trial function meth- 
ods and variational principles. The solutions are found to be consistent with the 
dynamical evolutions. 
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Chapter 1 



Introduction 



In our everyday lives, we typically think of bubbles as objects that separate two 
different phases of matter from one another (like the bubbles you see in boiling 
water that separate steam from liquid water). Within this work, however, a bubble 
is something a little more abstract: it is any scalar field configuration that interpo- 
lates between the minima of a double-well potential. These scalar field bubbles are 
still closely analogous to boiling water, where the minima of the double-well poten- 
tial act like the two different states of matter. Although an emphasis is placed on 
cosmological bubbles, the model can be applied to other branches of physics model- 
ing phenomena ranging from the polarization states in ferromagnets to topological 
defects in superfluids (actually, anything described by Landau-Ginzburg theory of 
phase transitions!). This work concerns itself primarily with the collapse and colli- 
sion of bubbles and to a special type of collapsed bubble known as an oscillon. 

The dynamics of these bubbles are governed by the classical fiatspace non- 
linear Klein-Gordon (nlKG) equation. There is a long history in mathematics and 



physics of trying to find new (non-trivial) solutions to nonlinear wave equations. 
One such type of solution that is of interest to many is the soliton. The first sci- 
entific discussion of solitons was due to J. Scott Russell and was published in the 
Report of the British Association for the Advancement of Science, in 1845. The 
report describes the creation of a surface wave in a narrow shallow water channel 
following the abrupt stop of a boat. Russell followed the wave on horseback un- 
til "after a chase of one or two miles [he] lost it in the windings of the channel." 
Although he dubbed these localized nonlinear waves "Waves of Translation", they 
later came to be known as solitary waves or solitonsQ. 

In the last fifty years, interest in solitons has been revived by many math- 
ematicians and physicists. Most particle physicists, for example, used to believe 
that for there to be bound (particle-like) states in a relativistic field theory, quan- 
tum theory had to be introduced. However, this is not the case, since solitons are 
stable bound states of a nonlinear field theory that, heuristically, exist through the 
balance between a nonlinear attraction and a tendency to disperse. Few scientists 
believe that quantum mechanics and field theory will ever be replaced entirely by 
solitonic interactions within classical field theory, but much research has gone into 
understanding how the existence of classical solitons implies the existence of a cor- 
responding quantum solution pH]. Although oscillons eventually disperse and are 
therefore not solitons (which are stable bound states), they do remain localized for 
large times and can pass through one another just like classical solitons. 

The background theory and history of the nlKG equation is the subject of 



^ Although some people |pO[| consider a soliton to be "any spatially confined and nondispersive 
solution to a classical field theory" , many others would call such a solution a solitary wave, reserving 
the term soliton to further include the ability for two such solutions to pass through one another 
with only a phase shift or time lag [p8|. 



Chapter 2. After presenting the nlKG equation, typical bubble initial data and the 
basics of bubble dynamics are discussed. Previous investigations of kink/antikink 
soliton interactions within the (1+1) dimensional A(/>^ model are then presented, and 
the chapter concludes with an introduction to oscillons (their behavior and history). 
Chapter 3 discusses the numerical techniques used throughout the thesis. 
Since finite difference methods are used extensively in the dynamic simulations a 
brief background is included. A section on dissipation is also included since the 
incorporation of dissipation is integral to the success of the numerical methods 
employed. The reader is then introduced to a new coordinate system that is used to 
solve nonlinear wave equations on a (relatively small) static lattice in one and two 
spatial dimensions. The chapter concludes with a brief motivation for presenting 



the coordinates in the "3+1" or ADM form, |[,[]16|. 

Chapter 4 contains the discussion of spherically symmetric oscillons. The fi- 
nite difference equations used are introduced and the testing of the code is discussed. 
A new resonant solution within the symmetric double-well potential (SDWP) model 
is presented; the mode structure of the solution is analyzed and a time scaling law 
is shown to be present for the critical (resonant) solution. Lastly, the threshold of 
expanding bubble formation is explored within the asymmetric double- well potential 
model and another (different) time scaling law is also shown to exist. 

Chapter 5 discusses axisymmetric evolution of the nlKG equation in the 
context of bubble collisions. The finite difference equations used are introduced and 
the testing of the code is discussed. The generation and testing of initial data is 
presented; the initial data is constructed from boosting two rest-frame oscillons (like 
those of chapter 4) at each other. Parameter space surveys are conducted and the 



threshold of expanding bubble formation is found to exhibit a time scaling law. 

Trial functions and variational approaches to finding critical oscillon solutions 
are the main ideas discussed in chapter 6. A set of generic ordinary differential equa- 
tions for critical non-radiative oscillon solutions used in chapter 4 are rederived using 
trial function methods. With more constrained ansatz (gaussian and hyperbolic se- 
cant functions) the same approach is used to directly obtain (ie. solving algebraic 
not differential equations) a few of the basic attributes of oscillons. 

Finally, chapter 7 concludes the thesis with a summary of what was accom- 
plished in this work. An appendix is also included that discusses some basic oscillon 
attributes (with units!): size, shape, lifetimes, and whether or not they are expected 
to form black holes. 

1.1 Notations, Conventions, and Abbreviations 

Unfortunately, this work mixes a few conventions from different branches of physics. 
The metric signature used, for example, is the one typical to (modern) general 
relativity, ( — h + + ). However, the units used throughout this thesis are Plankian, 
those commonly used by high-energy physicists where h = c = 1 (not those of the 
typical relativist where G = c= 11). If the lengths, lifetimes, and masses of oscillons 
considered throughout this thesis are left in terms of the dimensionful-scale in the 
model, m, with dimensions [i^]^^, the h is not needed. However, if one asks for the 
mass of a particular oscillon (particularly an early universe oscillon), one needs to 
include appropriate factors of h (or other appropriate dimensionful constant). With 
luck, any confusion that arises regarding how to reinsert units used can be dispelled 
by studying the examples in Appendix A. 



Also, the term "critical bubble" is used throughout the cosmology community 
to refer to a bubble (in a model with an asymmetric double-well potential) whose 
radius is at, or above, the threshold for expanding bubble formation. In other words, 
the term refers to a bubble large enough that the volume energy driving the field to 
the true vacuum is greater than the surface tension trying to collapse the bubble. 
Such bubbles will always expand and contribute to a phase transition. However, 
we choose to avoid the use of the word "critical" in such a generic way as it has 
quite a specific meaning in the study of critical phenomena, where it describes a 
solution that lies exactly on the threshold of the phase transition being considered. 
Therefore, instead of using the term "critical bubble", we tend to use "expanding 
bubble". 

The models discussed here are often referred to as (n+1) dimensional, (n+1) 
refers to a system with n spatial dimensions and 1 time-like dimension. The sym- 
metry of the model will also be included where possible to help distinguish between 
models with the same dimensionality but different symmetries, eg. (1+1) plane- 
symmetric and (1+1) spherically symmetric. 

Although defined throughout the thesis, we also note here the following fre- 
quently used abbreviations: 

• ADWP: asymmetric double- well potential 

• SDWP: symmetric double- well potential 

• nlKG: nonlinear Klein-Gordon 

• KG: Klein-Gordon 

• MIB: Monotonically Increasingly Boosted 



FDA: Finite Difference Approximation 
PDE: Partial Differential Equation. 
CN: Crank Nicholson 



Chapter 2 



Theory and Background 



This chapter presents a brief background of the history of the model to be studied 
throughout this thesis, the nonhnear Klein-Gordon model with double-well poten- 
tials. We focus on the work by Campbell et al [^, in the 1970's on the res- 
onant structure found in the (1+1) dimensional (plane-symmetric) kink-antikink 
scattering with the symmetric double-well potential, since many attributes of the 
plane-symmetric model are also found in the spherically symmetric case. Although 
interesting to this author for their non-linear (mathematical) behavior alone, brief 
descriptions of possible physical applications are included throughout. 

2.1 The Nonlinear Klein-Gordon Equation and Bubbles 

Put simply, this thesis is devoted to studying a special type of solution of the 
Nonlinear Klein- Gordon Equation, 

dV 



\m 




W) 



^ wan J. 

Figure 2.1: The symmetric double well potential, Vs{(l)), as a function of (j> (left). The 
SDWP has two global minima (degenerate vacua) at (a) and (c) and an unstable local 
maximum at (b). The field configuration shown (right) is a kink profile that interpolates 
between the two vacua. For r <C rmaii the field is at the vacuum point (c) on the potential, 
while for r 3> r^aii the field is at the vacuum point (a) on the potential; where constant 
and in the vacuum, the field has no energy. For r « r^uaii, the field interpolates between 
the vacua and must leave the vacuum; the wall therefore has both potential energy and 
gradient energy. In spherical symmetry the bubble will always collapse since there is only 
surface tension from the wall and no volume energy within the bubble (due to the degenerate 
vacua) . 

where cp = cl){f, t) and the potential is of the form V{(l)) = —m (f) + acj) + Xcf) , 
for m, a, and A constant. The potential and some sample initial data for a = 
(the symmetric double well) and for a 7^ (the asymmetric double well) can 
be seen in figures 2.1 and 2.2, respectively. The first use of the A(/>^ theory to 



discuss phase transitions is usually credited to Landau and Ginzburg, ||3g]. The 
theory is widely applicable and has been used to describe many types of phenomena, 



ranging from phase transitions in uniaxial ferrolectrics [46|, to phase transitions and 
defect dynamics within superconductors and superfluids [q^, to cosmological phase 
transitions resulting from the spontaneous breakdown of gauge symmetries ([pi|], 



], 13, 11, [||] m, and I3], for general treatment; and ||, H, H, Q, and 
], for oscillon related studies). Although in this thesis the particle physics and 
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Figure 2.2: The asymmetric double well potential, V{(j)), as a function of (f) (left). The 
ADWP has two local minima (non-degenerate vacua) at (d) and (f) and an unstable local 
maximum at (e). Vacuum point (d) is referred to as the false vacuum since it is not the 
global minimum and is unstable to fluctuations. Vacuum point (f) is referred to as the true 
vacuum since it is the global minimum and is stable to fluctuations. The field configuration 
shown (right) is a kink profile that interpolates between the two vacua. For r <C r^aii the 
field is in the true vacuum, while for r 3> r^jaH the field is in the false vacuum; where constant 
and in the vacuum, the field has no energy. For r « r^ijaU , the field interpolates between the 
vacua and must leave the vacuum; the wall therefore has both potential energy and gradient 
energy. In spherical symmetry their is competition between the surface tension trying to 
make the wall collapse and the volume energy (from the more energetically favorable true 
vacuum) trying to make the wall expand. The bubble's fate depends on the balance between 
the two forces. 



cosmology vernacular is mostly used, the results are applicable to any scalar field 

model described by the non-linear Klein-Gordon equation with a SDWP or ADWP. 

When describing a particle theory, a local minimum describes a vacuum 

state; the global minimum is the true vacuum while all other minima are false vacua. 



For the SDWP in figure 2.1, the sample initial data interpolates between the two 
degenerate vacuum states. For r > r^aii the field is in the vacuum state described 
by point a on the potential, V{(j)), while for r < r^aii the field is in the vacuum state 
described by point c. At and around r = Vi^aii, the field interpolates between the 
two vacuum states and is called a domain wall. Almost all of the energy of the field 
is concentrated in the wall. There is no potential energy or gradient energy in the 
field away from r^jaii, since it is in the vacuum; whereas in the wall the field leaves 
the vacuum, so it has both potential energy and energy due to the gradients of the 
field. If r is a Cartesian-like coordinate and the field is plane-symmetric then the 
wall is a planar domain wall. However, if r is a radial (spherical) coordinate, the 
wall has the shape of a spherical shell, and is often referred to as a bubble wall. The 
term bubble is in analogy to bubbles in fluids that are created during a change of 
phase, like gas bubbles forming in liquids, only instead of separating different states 
of matter, these bubbles separate different vacuum states of a particle theory. 

Although the term bubble is used for both SDWP and ADWP alike, the anal- 
ogy actually is best suited to the ADWP (depicted in figure |2^) . In the ADWP, the 
bubble wall separates the volume of space in the false vacuum from the volume of 
space in the more energetically favorable true vacuum. This is exactly like a super- 
heated liquid undergoing a phase transition to a gas, where the false vacuum is the 
liquid and the true vacuum is the gas. In the liquid-gas transition thermodynamic 



10 



fluctuations cause gas bubbles to form in the liquid; the fate of the bubbles is deter- 
mined by the competition between the surface tension in the wall and the volume 
energy within the bubble wall. For a large enough fluctuation it will be energetically 
favorable for the bubble wall to continue to expand, thus filling up space with the 
more energetically favorable gas. This is exactly what happens with vacuum state 
phase transitions, except instead of thermodynamic fluctuations, the bubbles are 
usually nucleated by quantum fluctuations. 

2.2 The (1+1) plane-symmetric (0^ — l) model 

The Klein-Gordon equation in (1+1) dimensions with the SDWP 

(where (/> = (j){x,t)) has been studied extensively over the last century. One of the 
most exciting discoveries was that the model supports solitary waves, stable localized 



solutions to non-linear fleld equations. For equation 2.2 these solitary waves take 
the form 

^{x, t)± = itanh f2i^^^\ (2.3) 

where 7 = l/\/l — v'^ and v is the velocity of the kink or antikink (plus or minus 
sign, respectively) . Much of the research of this model was motivated by the attempt 
to understand particle physics (meson) scattering experiments. In this context, 
Campbell et al.^^, showed that there is a resonant structure within the kink- 
antikink parameter space. The possible outcomes (varying depending upon the 
initial velocity of the kinks, ±Vi) were reflection, annihilation, or collapse to a long- 
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lived but unstable bound-state]^. A resonance between the translational mode of 
the colliding kinks and their individual internal shape mode vibrations gives rise 
to a very intricate structure for the endstate as a function of initial velocity. This 
structure was later shown (in the context of cosmological domain wall collisions) 
to be fractal in nature, ^. The fractal structure was found within "n-bounce 
windows", where the kinks collide, reflect off one another, and then collide again 
after radiating away energy (the process repeats itself n times). 

2.3 What is an Oscillon? 

The term oscillon has a few different meanings depending on context, but here it 
refers to a time-dependent spherically-symmetric coherent localized solution to a 
non-linear field equation (SDWP or ADWP) which is unstable but long-lived com- 
pared to the typical time-scale involved in the problem. Oscillons, originally called 
pulsons, were first studied numerically with the SDWP in 1977 by Bogolyubskii and 
Makhankov. They showed that for a wide range of initial data, the behavior of a 
collapsing bubble was characterized by three stages: 

• Collapse: The bubble wall collapses toward r = and the field oscillates 
irregularly while radiating a large amount of its initial energy. 

• Pseudo-stable Oscillations: The oscillon settles into a state where it remains 
localized (the location of the bubble wall is bounded) and the field oscillates 
about the stable vacuum, radiating very little energy. The term pseudo-stable 
is used because although the oscillons are unstable, they can last for thousands 



^Transmission is prohibited due to lack of energy conservation 
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of times longer than the time predicted by Hnear analysis^ (which also is 
roughly the oscillon's period of oscillation) . 

Dispersal: Eventually, an unstable shape-mode triggers the dispersal of the 
oscillon and the system is left in the original (r — > oo) vacuum state. 



Nearly twenty years later, Copeland et al. [24|, performed a much more thorough 
(and computationally rigorous) investigation of oscillons. Oscillons were shown to 
be extremely long-lived for a wide range of parameter values for both the SDWP 
and the ADWP (with varying degrees of asymmetry). The perturbative methods 
used provided an explanation for two properties of oscillons: A) the existence of 
a minimum radius for oscillon formation created from static bubble collapse, and 
B) the need for the initial energy of the field configuration to be above a certain 
threshold. 



Although Copeland, Gleiser, and Miiller |24|, were the first to really dissect 
oscillons and to start exploring how and why they behave the way they do, their in- 
vestigations (as good ones do) raised as many (or more) questions as they answered! 
In particular, Copeland et al. did not explore the fine structure of the parameter 



space as [10| and [|| did for the (1+1) plane-symmetric case, nor did they explore 
what effect non-spherical excitations play on the stability of oscillons. This work 
explores these two points and others that arose throughout the process. 



Although the true extent of their longevity was not shown convincingly until 1995, 
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Chapter 3 



Numerical Analysis and MIB 
Coordinates 



This chapter reviews the basic numerical methods (both new and old) employed 
throughout our research and motivates some of the choices made in the notation 
and form of the equations used. Although a complete description of finite difference 
equations, stability analysis, and dissipation are well beyond the scope of this thesis, 
a basic explanation of one and two dimensional finite difference and dissipation 
operators is provided. Finally, a new technique which solves problems associated 
with solving massive field equations on a lattice is introduced here. Implementation 
of this technique is then detailed in chapters ^ and ^. 

3.1 Finite Differences: Definitions and Notation 

Since the dawn of calculus in the 17th century, there has been a desire to be able 
to solve differential (and partial differential) equations. However, for most of the 
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time since Newton and Leibniz, the majority of mathematicians and physicists have 
been Umited to solving differential equations in closed-form or with various types of 
perturbative methods. These barriers have in large part collapsed since the creation 
of the computer in the latter part of the 20th century. Faster and faster computers, 
coupled with new numerical methods, continue to allow people to solve equations 
that before were far out of their reach. There are a great many numerical methods 
that have been developed for the solution of partial differential equations (PDEs), 
including finite differences, finite elements, spectral methods, and more. A brief 
explanation of finite differences is presented here, focusing on aspects most relevant 
to this research. (The following subsections are largely based on the class and lecture 
notes by Choptuik, 0, ||). 

3.1.1 Discretization 

For problems such as those considered here, finite differencing provides a very natural 
and straightforward route to the approximate solution of time-dependent partial- 
differential equations. A finite difference approximation (FDA) to a partial differ- 
ential equation (PDF) is obtained by replacing a continuous differential system of 
equations with a discrete system of approximate equations. The spacetime domain 
is represented by a discrete number of points on a static uniform mesh that are la- 
beled by Xfc and t", for integer k, n. The discretization scale is set by /i = Ax = AAt, 
and the error associated with the approximation should go to zero as h goes to zero. 
For a continuum differential system, 

Lu = f (3.1) 
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where L is a differential operator, u is the continuous solution, and / is the contin- 
uous source function, we use 

Lu = f (3.2) 

to denote its finite difference approximation, where L is a discrete difference opera- 
tor, u is the discrete solution, and / is the discrete source function. Throughout this 
chapter, ' will denote a discrete quantity; for clarity this notation will be dropped 
in subsequent chapters and discrete quantities will be recognized by their space and 
time indexes, u^. 

3.1.2 Residual 

It is useful to rewrite the FDA above as 

Lit- f = 0; (3.3) 

for a fully explicit scheme, this equation can be solved exactly. However, for iterative 
schemes, the right hand side of ( |3.3| ) will actually have a non-zero value that is 
representative of how well u solves the system. This leads to the definition of a 
residual 

f = Lu- /, (3.4) 

where u is the "instantaneous approximation" of n, i.e. that which converges to u 
in the limit of infinite iteration. Thus, f gives a measure of how well u satisfies the 
FDA, and an iterative scheme is said to be convergent if the residual is driven to 
zero in the limit. 
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3.1.3 Truncation and Solution Errors 

The truncation error of a finite difference approximation is defined to be 

f = Lu-f (3.5) 

where we note that the discrete operator L acts on the continuum solution u. f 
often cannot be obtained exactly since the solution u is not usually known. 
The solution error is defined to be 

e = u — u (3.6) 

and is the direct measure of how different the approximate solution u is to the 
continuum solution u. Although u must be known exactly to know the exact solution 
error, an approximate solution error (the solution error to leading order) can usually 
be obtained numerically (see 3.1.6| below). 



3.1.4 Consistency and Order of an FDA 

A difference scheme with discretization scale /i is a consistent representation of the 
continuous system if the truncation error goes to zero as the discretization scale 
goes to zero. Furthermore, the difference scheme is said to be p-th order accurate if 

limf = 0(/i^) (3.7) 

for an integer p. All of the schemes in this thesis are second-order, so hereafter it is 
assumed that p = 2. 
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3.1.5 Richardson Expandability 

For a centered difference sclieme L, witli discretization scale h, the solution u is 
related to the discrete approximation u by 

u = u + 62/1^ + 64/1*^ H (3.8) 

where the Cj are h-independent functions. A non-centered scheme will also have odd 
terms 63/1^, 65/1^, etc. Equation |3.8| is referred to as a Richardson expansion. 



3.1.6 Convergence 

A finite difference approximation is said to be convergent if and only if 

f^O as h^O. (3.9) 

Showing convergence is of prime importance in numerical analysis as it is a statement 
that the solution obtained numerically really approaches the continuum solution 
as the discretization goes to zero. A useful formula used throughout this thesis 
describes what is called (by Choptuik) the convergence factor, 

Cf ^ ^l^LZ^, (3.10) 

U2h - Uh 

where Uih, U2h-, and Uh, are the solutions with discretization scales 4/i, 2/i, and 
/i, respectively. Using the expansion ( |3.8| ) it can be shown that for second-order 
approximations C/ ^ 4 as h^ 0. 

3.1.7 Difference Operators 



Tables 3.1 and 3.2 contain the one and two dimensional second order difference 
operators to be used later. They can be derived by Taylor expanding the solution 
using the discretization scale as the expansion parameter. 
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Operator Definition Expansion 

A/< (-3< + 4<+i-<+2)/(2Ar) dru\l + 0{h^) 



A^< 


(3<-4<_i+<_2)/(2Ar) 


9rw " + O (/i2) 


A,< 


«+i-<-i)/(2Ar) 


a^M " + (/i2) 


A^3< 


«+i - <_i) / (rf+i - rU) 


dru'l+O {h^) 


A,< 


{u^' - U-) /At 


dtu ;+^ + O (/i^) 


Mr""? 


-edi.[6< + <_2 + <+2- 
4«-i+<+i)]/(16At) 


(Ar)^ d^u " + O (/i^) 


^r^r 


«+i+<)/2 


<+Uo(/.2) 



Table 3.1: One dimensional two-level FDA operators, h = Ar = A ^At. 

3.2 Dissipation and Stability 

Up to this point, there have been no comments made regarding the stability of 
numerically evolved solutions using finite differences. Much to the dismay of com- 
putational physicists, many of the difference schemes one uses to approximate solu- 
tions to partial differential equations do not work because they lead to unphysical 
growing modes. To understand this effect (and how to correct it) consider the 

discrete analog to the Fourier decomposition of the continuous function (j){k, t) = 

1 f"" 

— / e~* ^ (j){x,t)dx (for k continuous), 
\/2tt J—oo 

-. oo 

^""(0 = ^^ E ^-^^'"C. (3-11) 

V ^vr m=-oo 

for — TT < ^ < TT and discrete m, |64], |^]. A mode analysis can be performed 



by inserting the decomposition (3.11) into the difference equation and looking at 



the behavior of each Fourier mode over time. For linear differential equations, 
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Operator Definition Expansion 

^>l, (-3<,+4<+i.,-<+2,,)/(2Ai?) dRu\'l^+0{h^) 

AK. (3<. - 4<-i,, + <-2.,) / (2Ai?) afl^l"^. + O (h^) 

^rKj «+i,, - <-i,,) / (2Ai?) dRu\l^+0{h^) 



A^mJ^j 


(-3<^.+4z.^^.+i-<^.+2)/(2Az) 


9.K "^. + (h^) 


2 *J 


(3<,-4<,_i+<,_2)/(2Az) 


d.ul+0{h-) 


A^w"j 


«,-i-<,+i)/(2Az) 


d.u'" .+0 (h^) 


^*<. 


«;'-<,)/Ai 


dtu ;;+^ + o {h^) 


M^'^<, 


-edis[6<j- + <_2j + <+2,j- 
4 «-i,, +<+!,,)]/ (16Ai) 


{ARfd^^ul + 0{h^) 


f^f'Kj 


-edis[6<j- + u^j-_2 + <j+2- 
4«,-i+<,+i)]/(16At) 


{Azfdtul^+0{h') 


,,ave n 


(<r+<.)/2 


u "+^ + O (Af2) 



Table 3.2: Two dimensional two-level FDA operators, h — AR — Az — A ^At. 
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it is straightforward to write the sokition to the difference equation in the form 
u^~^^{k) = p (^) u^{k), where p is a complex function of (,. It is clear that if |/9p > 1 
the mode grows over time, if |pp = 1 the mode remains constant, and if |pp < 1 
the mode decays over time; |/op is referred to here as an amplification factor. A 
difference equation is dissipative if no modes grow with time and at least one mode 
decays, while it is nondissipative if the modes neither decay nor grow, or unstable if 
any of the modes grow with time. 

One method commonly used to make unstable schemes stable, or to make 
nondissipative and dissipative schemes more stable, is to add dissipation. Dissipa- 
tion can be added in a variety of ways, but we add it to our difference scheme by 
incorporating higher order spatial derivatives multiplied by the grid spacing to some 
power, (Ax)". Dissipation usually lowers the amplification factor for most modes, 
and typically affects higher wavenumbers more dramatically. Put another way, the 
goal is to dampen high frequency modes while maintaining the order of the original 
difference scheme. 

Using the dissipation operators in conjunction with second order CN deriva- 
tive operators (see table ^111 ), the linear advection equation {dt — d^) (j) = 0, can be 
approximated by the following FDA: 

dt 4dx Mx 

+t| {K + ^i+2 + uU - 4 {u]^2 + uU) } • (3-12) 

This scheme has an amplification factor |/0^|(C), where 

1 + i^ sin(0 - ^ (3 + cos(2e) - 4 cos(e)) 

piO = ^ • (3.13) 

l-i-sin(e) 
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Plotting IpP as a function of wavenumberQ, figure 3.1 shows that the CN scheme is 
nondissipative for e = and stable for e < 1. 

Unfortunately |pp(^) is not so easy to compute for general difference schemes. 
In particular, when the equation is nonlinear the stability of a difference scheme 
cannot be easily determined (in closed form), even if the derivative operators used 
are known to be stable in the linear case. We therefore take an empirical approach 
and include dissipation as needed in order to make the scheme stable. In fact, for 
evolution equations such as those studied here, the incorporation of dissipation is 
often essential to the construction of stable schemes. 

3.3 Geometry 

This work is certainly not the first to explore numerical solutions to massive scalar 
field equations. The literature is full of work on both massive {m?(jp') and nonlinear 
(e0^, A0^, sin((?i>), etc.) potentials in one, two, and even three dimensions, both cou- 
pled to gravity and in flatspace [^2|, [S3|, |1C]. One of the well-known problems with 



solving the nonlinear or massive Klein-Gordon (KG) equation (even in flatspace) is 
that there is no closed-form out-going boundary condition. The massive KG equa- 
tion is simply (□ — ml 4){x, t) = 0, which has a dispersion relation w^ = /c^ + m^ . 
Therefore, the velocity of the outgoing radiation cannot be uniquely determined 
and a satisfactory outgoing wave condition cannot be applied. While some of the 
radiation that reaches the outer boundary will leave the computational domain, sig- 
nificant amounts will also be reflected back and can contaminate the solution. If 
the phenomena being studied is short-lived, the computational domain can be made 



^Again, this is true only for the hnear advection equation. 
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Figure 3.1: Amplification Factors for CN schemes with and without dissipation, where e 
is the parameter measuring the amount of dissipation, e = corresponds to no dissipation 
while e = 1 completely quenches the modes at the Nyquist limit, ^ = tt. The CN scheme is 
marginally stable with no dissipation and stable for e > 0. 
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large enough so that radiation that does get reflected off the outer boundary will 
not have time to reach the region of interest. However, for long-lived phenomena 
(like oscillons or some boson stars) one must deal with the outgoing radiation more 
directly. 

Many previously used attempts consist of using an approximate out-going 
boundary condition and some sort of absorbing region near the outer edge of the 



grid. ||3[, 1 53 1, ||6^, and |6^] use a sponge filter, which imposes an outgoing radiation 
condition over a finite portion of the computational domain; this allows outgoing 
radiation to propagate while attempting to dampen ingoing radiation. However, due 
to the aforementioned unknown radiation velocity, this is done only approximately 
and can be susceptible to back-scatter effects. Recently Gleiser et al, [p9| ], have used 
a method referred to as adiabatic damping, where instead of focusing only on the 
outgoing radiation (as in the sponge filter methods) with a potential (roughly) of 
the form 7 (x) ((/) + (/)'], they use a term of the form 7 (x) cj) and dampen all scalar 
radiation. However, by adding explicit damping terms to the equations of motion 
(as opposed to higher order dissipation added to the difference equations) neither 
of the resulting difference schemes actually reduce to the true equations of motion 
in the limit that the grid spacing goes to zero. 

This section introduces a new geometric technique that effectively absorbs 
outgoing radiation, has difference equations that reduce to the differential equations 
in the continuum limit, and that is natural and straightforward to implement in both 
spherical and axial symmetry. The method employed has two parts, the transforma- 
tion of coordinates and the incorporation of dissipation into the numerical scheme. 
The coordinate system used leaves the interior of the grid alone, while transform- 
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ing the coordinates of the exterior points to be moving outward at approximately 
the speed of light relative to the interior or original rest frame; the coordinates 
are monotonically increasingly boosted (MIB). Characteristic analysis of the wave 
equation (in MIB coordinates) shows that both ingoing and outgoing characteristic 
velocities approach zero in a region near the outer edge of the grid [pC| ] . As the field 
slows down it becomes compressed; since the dissipation becomes stronger with in- 
creasing wavenumber, the field is quenched. This is shown to occur in a stable and 
non-reflective manner in (1+1) spherical symmetry and (2+1) axisymmetry. 

3.3.1 Radial MIB Coordinates 

In spherical symmetry, the outgoing radiation is frozen-out by introducing a new 
radial coordinate that smoothly interpolates between the standard polar radial co- 
ordinate f on Minkowski space 

ds'^ = -dp + df^ + f^dO^, (3.14) 

(where dCl'^ = dO'^ + sin 9 dcp^), and an outgoing null coordinate. We define 

i=t, f = r + f{r)t, and n = n (3.15) 

where /(r) is a monotonically increasing function that interpolates between and 
approximately 1 at some characteristic cutoff, Tc, 

{0 for r ^ Tc 
\ . (3.16) 

~ 1 for r S> Tc J 

In general, these coordinates are not good coordinates everywhere. However, if /(r) 
is monotonically increasing, the determinant of the Jacobian of the transformation is 
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non-zero for all t > — niax|/'(r)| and the coordinate transformation is one-to-one. 
Although a coordinate singularity inevitably forms as t approaches past timelike 
infinity, this has no effect on the forward evolution of initial data specified at f = 0. 
We must also demand that /(O) = to maintain the condition for elementary flatness 
at the origin. This coordinate choice takes the metric to 



ds^ = {-l + f{r))dt'^ + 2f{r){l + f'{r)t)dtdr 
+ (1 + f'{r)tf dr^ + (r + f{r)tf dft^ 

or in a more familiar (3+1) form ( p] , p^ , fig] ) to 



(3.17) 



ds'^ = (-a^ + a^/?^) dt'^ + 2a'^j3dtdr + a'^dr'^ + r'^b'^dn'^, (3.18) 



where 



air,t) = l + f'{r)t b{r,t) = 1 + /(r)- 

r 

«{M) . 1 /3(M) ^ :f^^ 



(3.19) 



Assuming the following specific form for / (see figure 3.2): 



/(r) = [1 + tanh ((r - r,) /5)] /2 + e(re, 6) , (3.20) 

where e{rc,5) = — [1 + tanh ((— Tc) /5)] /2, we can obtain the actual metric vari- 
ables, the characteristic velocities, and the conformal structure of the new hyper- 
surfaces. 

The characteristic analysis of the Klein-Gordon equation with metric ( 3.18| ) 
yields characterstics 

A± = -/3±-, (3.21) 

a 

where A+ and A_ are the outgoing and ingoing characteristics, respectively [^], [p^j. 

The MIB system behaves like the old {i,f) coordinates for r <C re, but the outgoing 
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Figure 3.2: Interpolating function, f{r'), for radial MIB coordinates, where r' is in units 
where re = 1. S is taken to he S ^ 6/rc ~ 0.0893 (corresponding to the system to be used 
in chapter 4). Wherever / = 0, the new radial coordinate is left unchanged with respect 
to the old coordinate, r. Where / > the new radial coordinate is being "shifted" with 
respect to r with velocity (3 = //(I + ft). 
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Figure 3.3: Plot of the characteristic velocities as a function of r' and t', where r' and t' 
are radial MIB coordinates in units where Vc is set to unity, and A+ and A_ are the outgoing 
and ingoing characteristics, respectively. 8 is taken to be (5 ^ i5/rc ~ 0.0893 (corresponding 
to the system to be used in chapter 4). Characteristic velocities are plotted for times t' = 0, 
1, 10, and 100 (i' = 100 is larger than the lifetime of the longest lived solution studied in 
this work). Both characteristic velocities approach zero around r' = 1 as 1/t'. 
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Figure 3.4: Plot of characteristic velocites, A±(r', 100), where r' and t' are radial MIB 
coordinates in units where re is set to unity, and A+ and A_ are the outgoing and ingoing 
characteristics, respectively. S is taken to he 5 ^ S/tc ~ 0.0893 (corresponding to the 
system to be used in chapter 4), and t' = 100 is larger than the lifetime of the longest lived 
solution discussed in this thesis. Both characteristic velocities approach zero around r' ~ I 
as 1/t'. 
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radiation gets frozen out in the r ~ Tc region. Figure 3.3 and equations ( p. 21 ) and 



( |3.19| ) show that around r ^ r^ both the ingoing and the outgoing characteristic 
velocities go to zero as t ^ oo (as the inverse power of t). It is this property that 



is responsible for the "freezing-out" of the outgoing radiation |2C]. We call these 
coordinates (t,r) monotonically increasingly boosted (MIB) radial coordinates. 

The conformal structure (figure ^]^) is obtained by applying equations (|3.15|) 
to the standard conformal compactification on Minkowski space, t it f = tan ( ^ ) 
(where here T and R are the axes in the conformal diagram, see [^] or |6g] ) , and then 
plotting curves of constant r and t. The constant-t hypersurfaces are everywhere 
spacelike, and all reach spatial infinity, i°. Although constant-r surfaces for r > r^ 
appear at first glance to be null (or as if ( ^1 "tips over" to become a null vector), 
closer examination (see insets of Fig. |3.5| ) reveals that they are indeed everywhere 
timelike and do not reach future null infinity, J'^ . This behaviour follows from the 
fact that the interpolating function, /(r), only reaches unity asymptotically, at i°. 

3.3.2 Axi-Symmetric (2-D) MIB Coordinates 

The (1-D) MIB coordinates introduced in the previous subsection also have a straight- 
forward implementation in two dimensions. The (2-D) problems discussed in this 
work all model (3-|-l) dimensional spacetimes with axial symmetry. (2-D) MIB co- 
ordinates are obtained by again starting from the Minkowski metric, but now in 
cylindrical coordinates 

ds^ = -dp + dR^ + R^de'^ + dz^. (3.22) 

These coordinates are modified just like the radial coordinate in spherical symmetry 



(section |3.3.1 ), except now both the axisymmetric coordinates R and z interpolate 
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Figure 3.5: Conformal diagram sliowing surfaces of constant r (dotted lines) and lines of 
constant t (dashed lines). Lines of constant t look exactly like the constant-t hypersurfaces 
of Minkowski space, whereas the lines of constant r behave much differently. For r > re, it 
appears as if the constant-r surfaces are null. This occurs since the coordinates are being 
shifted outwards at nearly the speed of light. However, as insets A, B, and C show, the 
constant-r lines do not become null (do not intersect future null infinity), and are everywhere 
timelike. This is actually easy to understand as the radial coordinate is never moved out at 
the speed of light (not even at spatial infinity, because of the e term in equation 3.20|). 
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between the old {R,z) coordinates and "null coordinates at i°": 

i=t, R = R + f{R)t, z = z + g{z)t, and 6 = 6. (3.23) 

Here, f{R) and g{z) are monotonically increasing functions: f{R) interpolates be- 
tween (approximately) and 1 at a characteristic cutoff, Re, g{z) interpolates (ap- 
proximately) between —1 and 1 at characteristic cutoffs, i^^cQi 

— 1 for z <^ —Zc 



for R<^ Re 

1 for ii > Re 



for —Zc < z < Zc 

1 for z ^ Zc 



(3.24) 



This coordinate choice takes the metric to 

ds"^ = [-1 + f{Rf+g{zf)dt^ + 2f{R){l + dRf{R)t)dtdR 
+2g{Z) (1 + d^g{z)t) dtdz + (1 + dRf{R)tf dR'^ 
+ {R + f{R)tf dO^ + (1 + d,g{z)tf dz^ 

or in a more familiar (3-|-l) form to 

ds2 = (-a2 + a2/3^^ + &2^^2^ dt^ + 2a'^f3^dtdR + 2b'^(3^dtdZ 



(3.25) 



WdR^ + R^dO'^ + h'^dz'^, 



(3.26) 



where 



a{R,t) 

f3''iR,t) 

a{R,t) 



l + dRf{R)t 

fJR) 
1 + dRf{R)t 

1 



b{z,t) = l + d,g{z)t 



f3'{z,t) 



R{R,t) 



1 + d,g{z)t 
R + f{R)t. 



(3.27) 



In accordance with ( |3.24| ), the interpolating functions, / and g, are taken to be 

f{R) = [l + tanh{{R-Re)/6R)]/2 + e{Re,6R) (3.28) 



9iz) 



tanh {{z - Zc) j^z) /2 + tanh ((z + z^) /d^) /2, 



(3.29) 



^We could have chosen to have to distinct cutoffs, {zc)^ and {zc) , but by performing the 
"experiment" around z = this simphfication works perfectly well. 
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where e(i?c, Sr) = - [1 + tanh {{-Re) /Sr)] /2, (see Figure U). 
The characteristics of the system are 

A^ = -I3^±l/a (3.30) 

A| = -p^±l/b, (3.31) 

where A^ and A^ are the characteristics in the R and z directions, respectively. As 
with the spherically symmetric case, the system behaves like the old {R, z, and t) 
coordinates for R < Re and —Zc < z < Zc, while around R k, R^ and z ~ i^c, the 
characteristic velocities, A^ — > and A^ — > 0, as t ^ oo (as the inverse power of t). 
This again has the effect of freezing-out the outgoing radiation. 

3.4 Dissipation &; MIB Coordinates Working Together 

The primary reason dissipation was introduced and stressed in this chapter is be- 
cause it is integral to the stability of the MIB scheme. Since MIB coordinates cause 
all the outgoing radiation to accumulate in a very small region of the grid (creating 
large gradients in the fields), when MIB coordinates are used without dissipation 
the scheme is always unstable. However, with dissipation the field is quenched in a 
stable and non-reflective manner. This can happen because the amplification factor 



for the Crank-Nicholson scheme (analagous to that seen in figure 4.2) is significantly 
less than one for high frequency solution components. Therefore, the more the ra- 
diation gets spatialy compressed, the more it is dissipated. Details of this process 
are discussed in chapter ^ in the context of spherically collapse. 

Note that although the above explanation is self-consistent and explains how 
the system can be stable, it certainly provides no convincing argument that the 
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scheme shouldhe stable (or non-reflective for that matter). In fact, many techniques 
explored throughout the creation of the one discussed here were very unstable and 
reflected significant amounts of radiation! In addition to the sponge-filter approach 



(mentioned in section 3.3), there has been much research devoted to the study of 
outgoing radiation and outgoing boundary conditions. For the sake of completeness, 
we briefly review some of this research here. 

Cauchy-Characteristic Matching (CCM) is a way to solve Einstein's Equa- 
tions that matches typical Cauchy evolution on the interior of the computational 
domain to characteristic evolution on the exterior, [Q]. The original motivation was 
to measure the outgoing radiation that reaches null infinity, however, in the pro- 
cess the method eliminates the need for an outgoing radiation boundary condition. 
The technique was developed for massless radiation although some success has been 



made recently with the incorporation of matter H], |51]. CCM is reported to be 
more efficient than previously used techniques, but still requires separate evolution 
of two domains and calculations to match them together; the methods employed 
here involve a single Cauchy evolution. 

Null-cone evolution |^l| and other purely characteristic evolution techniques 



(ie. where the entire evolution is performed in characteristic coordinates, |40|, ^ 
have also been successful in numerical relativity and seemed like a very promising 
way to evolve the nlKG equation. Unfortunately, using characteristic coordinates 
with a truncated grid is still problematic due to the lack of a satisfactory outgoing 
boundary condition. However, unlike typical (r,t) coordinates, compactification is 
natural with null evolution and generally preserves the smoothness properties of ra- 
diation [B^. Furthermore, in compactified characteristic coordinates exact boundary 
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conditions can be set on the outer boundary of the grid, J^ . Numerical implementa- 
tions can usually be evolved stably and radiation can be measured at J^ . However, 
this tends to work well only for radiation that travels along null characteristics (ie. 
massless fields). When the fields are massive the field does not propagate along 
null characteristics (and therefore does not reach null infinity), and there is either 
leakage out to J'^ of massive field, or instabilities that arise from steep gradients 
due to compression of the outgoing radiation, [[7C| ]. 

Similar problems occur evolving massive radiation with the conformal com- 



pactification methods developed by Friedrich, |34], and implemented numerically by 



Hubner, |44] and Frauendiener, |28|, |29|, ]3^; the conformal compactification of the 
spacetime brings null infinity to the outermost grid point and still causes problems 
with massive radiation [^]. 
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Figure 3.6: Interpolating functions, f{R') and g{z'), for axisymmetric MIB coordinates, 
where R' and z' are in units where Re and Zc are equal to one. 6fj and S^ are taken to 
be approximately 0.1067 (corresponding to the system to be used in chapter 5). Wherever 
f — 0, the new radial coordinate, R, is left unchanged with respect to the old coordinate, 
R. For / > the new radial coordinate is being shifted outward with respect to R with 
velocity f3^ — //(I + daft). Wherever g{z') = 0, the new z coordinate is left unchanged 
with respect to the old coordinate z. For <? > the z coordinate is being shifted outward 
with respect to z with velocity f3^ = g/{l + dzgt); note however that the sign of g (and 
consequently /3^) is negative for z < as the coordinate is moving outward toward more 
negative z. 
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Figure 3.7: Plot of the characteristic velocities, A^, as a function of B! and i', where B! 
and t' are axisymmetric MIB coordinates where Re is set to unity, and A+ and A_ are the 
outgoing and ingoing characteristics, respectively, b is taken to be 8r -^ 5r/Rc ~ 0.1067 
(corresponding to the system to be used in chapter 5) . Characteristic velocities are plotted 
for times t' = 0, 1, 10, and 100, and where t' = 100 is larger than the lifetime of the longest 
of the longest lived solution studied in this work. Both characteristic velocities approach 
zero around r' = 1 as 1/t' . 
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Figure 3.8: Plot of the characteristic velocities 
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as a function of 



and t' where 



and t' are axisymmetric MIB coordinates where Zc is set to unity, and A+ and A_ are the 
outgoing and ingoing characteristics, respectively. 5z is taken to be 5z -^ 5z/ Zc ~ 0.1067 
(corresponding to the system to be used in chapter 5). Characteristic velocities are plotted 
for times t' = 0, 1, 10, and 100 (t' = 100 is larger than the lifetime of the longest lived 
solution studied in this work). Both characteristic velocities approach zero around z' — ±1 
as l/t'. 
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Figure 3.9: Plot of the characteristic velocities 
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as a function of 



and t' where 



and t' are axisymmetric MIB coordinates where Zc is set to unity, and A+ and A_ are the 
outgoing and ingoing characteristics, respectively. 5z is taken to be 5z -^ 5z/ Zc ~ 0.1067 
(corresponding to the system to be used in chapter 5). Characteristic velocities are plotted 
for times t' = 0, 1, 10, and 100 (t' = 100 is larger than the lifetime of the longest lived 
solution studied in this work). Both characteristic velocities approach zero around z' — ±1 
as l/t'. 
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Figure 3.10: Plot of characteristic velocites, \±{R' , 100), where R' and t' are axisymmetric 
MIB coordinates in units where Re is set to unity, and A+ and A_ are the outgoing and 
ingoing characteristics, respectively. 5 is taken to be (5 ^ 5/ Re ~ 0.1067 (corresponding to 
the system to be used in chapter 5), and t' — 100 is larger than the lifetime of the longest 
lived solution discussed in this thesis. Both characteristic velocities approach zero around 
r' = 1 as l/t'. 
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Figure 3.11: Plot of characteristic velocites, A±(z', 100), where z' and t' are axisymmetric 
MIB coordinates in units where Zc is set to unity, and A+ and A_ are the outgoing and 
ingoing characteristics, respectively. 5 is taken to be (5 ^ (5/zc ~ 0.1067 (corresponding to 
the system to be used in chapter 5), and t' — 100 is larger than the lifetime of the longest 
lived solution discussed in this thesis. Both characteristic velocities approach zero around 
r' = 1 as l/f. 
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Chapter 4 



Spherically Symmetric Oscillons 



4.1 The Klein-Gordon Equation in MIB Coordinates 
with SDWP 

In this chapter we introduce the notation, equations of motion, and discuss results 
from a computer code which uses spherically symmetric MIB coordinates. 

We are interested in massive scalar field theory described by the (1+1) spher- 
ically symmetric action 



5[0] = I d^x^\ (^-]^g^-y^<pv,4> - V{<p)) (4.1) 



where (/> = cj){r,t), V^(0) is a symmetric double well potentialF] (see figure |4.1| ), 
Vs{cl)) = — ((/) — l) , dfiu is the flatspace metric in spherically symmetric MIB co- 
ordinates, (p. 18), and g is the determinant of g^^- The equation of motion for the 



^This is identical to using V{(j)) = 7 I <^ ^ ) and introducing dimensionless variables r = 



fm, t = tm, and Tp — —q 
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> 




Figure 4.1: Symmetric double well potential, Vsicj)) = - 



-1 



The potential has two 



degenerate minima at (p — ±1 and an unstable local maximum at = 0. See figure |J| for 
detailed description of features and interpretation of kink profile (bubble) initial data. 
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action ( ^ ) is 



5m ^\9\9''''d. 



which with ( |3.18| ), (|3.19 ), and the definitions 



n(r,t) = -{dtct>-(idA), 



a 



$(r,t) = dr^, 



yields 



n 
6 



a 



^252 f»^ + pu]] - 2^n - aacp (c/.^ - 1 



a 



n + /3$ 



a 



n + /3$ 



(4.2) 

(4.3) 
(4.4) 

(4.5) 

(4.6) 
(4.7) 



where ' = dt and ' = dr- These equations are the famihar (3+1) form for the 
spherically symmetric Klein-Gordon field coupled to gravity. However, instead of 
having a truly dynamical geometry, here a{r, t), b{r, t), a{r, t), and /3(r, t) are known 
functions of (r, t), which result from the MIB coordinate transformation of flatspace. 



4.2 Finite Difference Equations 



Equations (^3,4.(:.4?7|) are solved using two-level second order (in both space and 
time) finite difference approximations on a static uniform grid with Nr grid points. 
The scale of discretization is set by Ar and At = AAr, where we fixed the Courant 
factor. A, to 0.5 as we changed the base discretization. 

Using the operators from Table ^.1| , drV = a, dr = nr^~^dr", and rb = f, the 
interior difference equations (3 < i < A^^ — 2) are 



Atn- 



3Att 



aAf3 [r' [-^ + I3U 
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-2f^J^U-aa<p{cp'-l)\ + ^f ^H^, 



(4.8) 



A,$^ = /itA, (^n + /3cDJ" + ^di-$«, (4.9) 



At</.^ = ;,, f-n + /3$l +/if^^(/.^ (4.10) 

and include dissipation. On the two grid points adjacent to the boundaries, i = 2 
and i = Nr — 1, the same equations are solved without dissipation (since the two 
points needed on either side are unavailable): 



AiH^ = 3fit 






-2fit { -n - aa(t> [4>^ - 1) \ , (4.11) 



a 



At^': = fitAr ^-U + /?<!> j , (4.12) 

AtC = fit(^Il + ^^y. (4.13) 

These equations are solved using an iterative scheme, which is iterated until the 
^2-norms of the solutions at t"^^^ converge to one part in 10^, where the ^2-norm is 
defined for a vector Xj to be 

/ -. N \ 1/2 

Since the functions a, b, b, a, /3, and f are explicitly known functions, when dis- 
cretized we simply evaluate the given function at Vi = (i — l)Ar and t" = (n — l)At. 
On the inner (r = 0) boundary we applied conditions necessary for regularity at the 
origin: 

$^ = (4.15) 

/ijAfn-^nj =0 (4.16) 
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where ( |4.15[ ) is a statement of regularity in the scalar field, and ( 4.16| ) results from 



equation ( [4.13 ) and the commutation of partial derivatives. The field (p is evolved 
using 

Ai</. = /Xi(^— + /3$j . (4.17) 

For the outer boundary, our choice of equations makes very little difference since 
the physical radial position (f) corresponding to the outermost point is moving 
out at nearly the speed of light, therefore none of the outgoing field ever reaches 
this gridpoint. Nevertheless, we employed the typical massless scalar field outgoing 
boundary condition for IT and $, and evolved <J3 with its equation of motion: 

/ rr\ " 

A^n- + f,t (^Af,n + - j . = 0, (4.18) 

At<^^ + fit l^/^t'^ + ^y = 0, (4.19) 

/all \" 

AtcP = f,t^— + p^j . (4.20) 

4.3 Testing the MIB Code 

One might think that freezing out the outgoing radiation while keeping a static 
uniform mesh would lead to a "bunching-up" of the outgoing radiation from the 
oscillating source which would cause a loss of resolution, numerical instabilities, and 
eventually crash the code. However, as already discussed, this turns out not to be 
the case; all outgoing radiation is stably "frozen-out" around r ~ Tc and the steep 
gradients that should form in this region are smoothed out and quenched by the 
dissipation (see figure |4.2| ). As the oscillations approach r th r^ their wavenum- 
ber approaches the Nyquist limit, ^ — > vr, and since the amplification factor for 
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these modes is significantly less than one they are rapidly quenched (analogously to 



the solution of the advection equation using the CN scheme shown in figure 3.1). 
Put simply, the coordinate system causes the oscillations to "bunch-up", while the 
higher-order dissipation quenches the resultant high-frequency modes. 

There is a loss of resolution and second order convergence directly around 
Tc, but this does not affect stability or convergence of the solution for r <^ re- 
Figure ( [4. 3D shows a convergence test for the field (/) for r < rc/2 over roughly six 
crossing times. Since we are solving equation ([4.2| ) in flatspace, it is very simple to 
monitor energy conservation. The spacetime admits a timelike Killing vector, t'^ , 
so we have a conserved current, J^ = f^Tf^y. We monitor the flux of J^ through a 
gaussian surface constructed from two adjacent spacelike hypersurfaces for r < r^ 
(with normals n^ = (±1,0,0,0)), and an "endcap" at r = rendcap (with normal 
n^ = (0, a^^,0, 0)). To obtain the the conserved energy at a time, t, the energy 
contained in the bubble, 

Ebnhhieit) = 4vr I rH^ I ^^, + V{(P) J dr, (4.21) 

^ ^ 

(where the integrand is evaluated at time t) is added to the total radiated energy, 

t 

r „ „T\<b 

ErUt) = A-K r%^—dt' (4.22) 

J Hi 



(where the integrand is evaluated at r = rendcap)- The sum, -Etotai = -E'bubbie + -E'rad, 

remains conserved to within a few tenths of a percentfl through almost two thousand 

oscillation periods {t ~ 8000m^^ or up to a quarter million iterations, see Fig. ( [4.4] )). 

Although monitoring energy conservation is a very important test, it says 

little about whether there is reflection off of the outer boundary or the r ~ Tc 

^A few hundredths of a percent if measured from after the initial radiative burst from the 
collapse. 
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55 40 



Figure 4.2: Fundamental field 0(r) in the freeze-out region (cutoff actually at Tc — 56) 
at t = 0, 57, 75, and 27(1 As the characteristic velocities of the radiation go to zero 



around Tc (as seen in figure 3.4), the wavelength of the radiation is shortened to the Nyquist 
limit on the lattice, 2Ar. The higher-order dissipation added to the system is responsible for 
quenching the field; figure 3.1 shows the amplification factor \p^ as a function of wavenumber 
is significantly less than one as the wavenumber approaches the Nyquist limit (the plot of 
|p(f)P is actually for the advection equation, but is qualitatively similar). 



48 




Time 



Figure 4.3: Convergence factor, C/ = \\(j)4h — (t>2h\\2 / \\4'2h — <f>h\\2, for the field (j) com- 
posed from the solution at three different discretizations (value of 4 indicates 2nd order 
convergence, as shown in section 3.i.6). 



49 



150 



100 - 



Si 
CD 




50 - 



2000 3000 
Time 



4000 



5000 



Figure 4.4: Plot of energy contained in oscillon (dashed green lines), energy radiated 
(dotted red lines), and total energy (solid blue line). The total energy of the system is 
a constant of motion and is numerically conserved to within a few tenths of a percent 
(hundredths of a percent if measured from after the initial radiative burst from the collapse). 
The energy contained within the oscillon drops rapidly during the initial radiative phase and 
plateaus around E w A'irn/X during the pseudo-stable regime of the oscillon. 
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region. To check this we compare output from the MIB code to two other types 
of numerical solutions. The first type of reference solution involves evolution of 
equation ( |4.2| ) in {r,i) coordinates, but on a grid so large that radiation never reaches 
the outer boundary (large-grid solutions). These solutions serve as our "ideal" 
reference solutions (for a given discretization), since we are guaranteed that they 
are free of contamination from reflections off the outer boundary. The second type of 
reference solution involves evolution on spatial grid comparable in extent to the MIB 
grid, but using an outgoing boundary condition (OBC solutions). Since we know 
these solutions do have error resulting from reflection off of the outer boundary, they 
demonstrate clearly what can go wrong when the solution becomes contaminated 



by reflected radiation. Treating the large-grid solution as ideal. Fig. (4^) compares 
typical log]^Q \\(j) — i?i>idoai||2 for the MIB and OBC solutions. There is a steep increase 
in the OBC solution error (three orders of magnitude) around t = 125, which is at 
roughly two crossing times. This implies that some radiation emitted from the initial 
collapse reached the outer boundary and reflected back into the region r < rendcap- 
There is no such behavior found in any MIB solutions. Lastly, for a more direct look 
at the field itself, we can see (j){0,t) for large-grid (triangles), MIB (solid curves). 



and OBC (dashed curves) solutions in Fig. (4^). Initially, both the MIB and OBC 
solutions agree with the large-grid solution extremely well, while after two crossing 
times the OBC solution starts to diverge. 

Since the MIB solution conserves energy to second order, converges quadrat- 
ically (in the domain of interest), and is equivalent to the large-grid (again, within 
second order error) , it is an acceptable means of solving equations ( [4.5| , [4.(]|j4^ ) while 
being a great deal more computationally efficient than other previously used tech- 
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Figure 4.5: Plot comparing the OBC (solid red) and MIB (dashed green) solutions to 
an "ideal" solution. The OBC solution is obtained using a massless outgoing boundary 
condition, the MIB solution is obtained by solving the system in spherical MIB coordinates, 
and the ideal solution is obtained by evolving the solution in standard (r,t) coordinates on 
a grid large enough to ensure no reflection off the outer boundary. The error estimates are 
obtained from the ^2-norm of the difference between the trial solutions (OBC or MIB) and 
the ideal solution, \\(f> — (/>idoai||2- Contamination of the OBC solution is observed at two 
crossing times, t w 120, where the error estimate increases over three orders of magnitude. 
The MIB solution error grows slowly and steadily as expected when solving a continuum 
equation with two different finite difference techniques. 
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Figure 4.6: 0(0, t) versus time for the large-grid solution (blue triangles), MIB 
solution (green solid curves), and OBC solution (red dashed curves). The solutions 
all agree before 2tcrossing) but the OBC solution begins to drift away from the ideal 
solution after 2tcrossing- The error in the OBC solution is due to radiation that is 
reflected off of the outer boundary (hence needing two crossing times to return to 
r = to contaminate the oscillon). All pictures span the same area, A(f) = 0.075 by 
At = 0.5. 
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niques. In fact, assuming that dynamical grid methods add more points hnearly 
with time, Nr{t) oc t, the total computational work, W oc tNr{t), grows as t^. With 
the MIB method, which uses a static uniform mesh, the computational demand 
grows linearly with the oscillon lifetime, W oc tN^ for A^^ constant. Thus, partic- 
ularly for long integration times, there is signifcant cost benefit in using the MIB 
system instead of a dynamically-growing-grid technique. 

4.4 The Resonant Structure of Oscillons 

{ID Critical Phenomena I) 



Copeland et al. [24|, showed quite clearly that oscillons formed for a wide range of 
initial bubble radii. By collapsing bubbles of many different radii, they even caught 
a glimpse of resonant structure (which in large part motivated this study) but they 
did not explore the parameter space in detail. With the efficiency of our new code, 
we are able to explore much more of the parameter space for a given amount of 
computational resources than with either large or dynamically growing grids. 

Following |2^ we use a gaussian profile for initial data where the field at the 
core and outer boundary values are set to the vacuum values ((/>c = 1 and (/>o = — 1 
respectively) and the field interpolates between them at a characteristic radius, rg: 

(t){r, 0) = 0o + {(pc - <Po) exp [-r^/rl) . (4.23) 

By keeping tpc and (po constant but varying tq, we have a one parameter family 



of initial data to explore. Figure 4.7 shows the behavior of oscillon lifetime as a 
function of tq. We discuss three main findings that are distinct from previous work: 
the existence of resonances and their time scaling properties, the mode structure of 
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Figure 4.7: Plot of Oscillon lifetime versus initial bubble radius for 2.0 < tq < 5.0. With 
the high resolution parameter space survey the 125 resonances become apparent. The survey 
was conducted at (relatively) coarse resolution to reveal the location of the resonances. Once 
a resonance was found, it was resolved cfhcicntly first by using a three-point routine that 
maximizes the lifetime as a function of initial bubble radius, and then by bisecting (to one 
part in 10^^) on the bifurcate modulation behavior discussed below (and seen in figure 4.9). 
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the resonant solutions, and the existence of oscillons outside the region 2 < tq < 5. 



4.4.1 Time Scaling 



A new feature is observed in the hfetime profile of collapsing bubbles; figure 4.7 
displays the appearance of 125 resonances. These resonances (also seen in Fig. 
[4.8| ) become visible only after carefully resolving the parameter space. Upon fine- 
tuning to one part in ~ 10^^ we noticed interesting bifurcate behavior about the 
resonances (figure |4.9| , top). The field oscillates with a period T ~ 4.6, so the 
individual oscillations cannot be seen in the plot, but it is the modulation that 
is of interest here|. The top figure shows the envelope of 0(0, t) on both sides of 
a resonance (dotted and solid curves). We see that the large period modulation 
that exists for all typical oscillons disappears late in the lifetime of the oscillon as 
tq is brought closer to a resonant value (which we define to be Tq). On one side 
of Tq the modulation returns before the oscillon disperses (solid curve), where on 
the other side, the modulation does not return and the the oscillon just disperses 
(dotted curve). This behavior does not manifest itself until tq is quite close to 
Tq. So in practice, we used a three point maximization (bracketing interval of 
~0.62) routine to get r^ close to Tq and then we bisected on the bifurcate behavior 
thereafter (bracketing interval of 0.5). This is much more efficient than a uniformly 
sampled parameter space; if the parameter space were sampled uniformly at the 
finest resolution used in figure 4.7, the survey would contain approximately 3 x 10^^ 
points, whereas we used only 7605. Although we can see from Fig. |4.9| that the 

modulation is directly linked to the resonant solution, it is not obvious why this 

^This would be r = 4.6m^^ in the original coordinate system. To recover the proper dimensions, 
lengths and times are multiplied by m~^ and energies by mX~^ . 
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Figure 4.8: Plot of oscillon lifetime versus initial bubble radius for 2.27 < riniuai ^ 
2.29. The three resonances shown occur at Tq ~ 2.2805, rl ~ 2.2838, and r^ ~ 
2.2876. Each resonance separates the parameter space into regions with n and n + 1 
modulations; the red points correspond to oscillons with no modulations, the green 
points to oscillons with one modulation, the blue to two modulations, and the violet 
to three modulations. 
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Figure 4.9: Top plot shows the envelope of (j){0, t) for rgiAro displaying bifurcate behavior 
around the Tq ~ 2.335 resonance (Aro ^ 10"^'*); the solid red curve is the envelope barely 
above resonance while green dotted line is the envelope barely below resonance. Bottom plot 
shows the energy radiated out of the gaussian surface containing the oscillon as a function 
of time. The increases in the energy radiated are synchronized with the modulation in the 
field. 
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is so. However, if we look at the relationship between the modulation in the field 
(top) to the power radiated by the oscillon (bottom), we see that they are clearly 
synchronized. 

This should really be no surprise for one familiar with the 1+1 KK scat- 
tering studied with the same model [10|. Campbell, et al., showed in a study of 
KK collisions that after the "prompt radiation" phase (the initial release of ra- 
diation upon collision) the remaining radiation emitted is from the decay of what 
they referred to as "shape" oscillations. The "shape modes" were driven by the 
contribution to the field "on top" of the K and K soliton solutions. Since the exact 
closed-form solution for the ideal non-radiative KK interaction is not known, initial 
data are only an approximation and the "leftover" field is responsible for exciting 
these shape modes. The energy stored in the shape modes slowly decays away as 
the kink and antikink interact and the solution eventually disperses. 

In our case, we believe the large period modulation is signaling the excitation 
of a similar "shape mode" on top of a periodic, non-radiative, localized, oscillating 
solution. On either side of a resonance in the rg parameter space, the solution is on 
the threshold of having one more shape mode oscillation. If this is the case, and we 
are "tuning in" this unstable shape mode, we might expect to see a scaling law arise. 
Figure [4.10| shows a plot of oscillon lifetime versus ln|ro — rg| (for the tq ~ 2.335 
resonance) and we can see quite clearly that there is a scaling law, T ~ 7 In |ro — TqI , 
for the lifetime of the solution on each side of the resonance. We denote 7hi for 
the scaling exponent on the r^ > Tq side, and 7L0 for the scaling exponent on the 
ro < Tg side. A plot of both scaling parameters as a function of rg is seen in Fig. 



4.11. We observe that for almost all of the resonances that 7HI ~ 7lo- 
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Figure 4.10: Plot of time scaling, T versus In |ro 
top (blue) line displays the scaling behavior for ro 



> r, 



for the ro « 2.335 resonance. The 
while the bottom (green) line for 



ro < rg. The exponents (measured by the slope) are both approximately equal to 7 = 30. 
Although each resonance has a scaling law, the exponents vary from resonance to resonance; 
a plot of the scaling exponent, 7, versus the critical initial bubble radius can be seen in figure 
Oil 
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Figure 4.11: Plot of critical exponents for each resonance. There are two values of 7 
for each resonance. The green values are the 7hi values measured on the tq > r^ side of 
the resonance while the blue values are the 710 values measured on the rg < Tq side of the 
resonance. The uncertainties are estimated from running the entire parameter space surveys 
at two resolutions and estimating the error, A7 = |7(Ar^=i449) — 7(Ar,,=io25)|- 
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4.4.2 Mode Structure 



If there does exist a periodic non-radiative solution to equation 4.2, we should be 
able to construct it by applying an ansatz of the form 



0(r,t) 



-1 + 



+ X! (f>nir) COS {nut) 



(4.24) 



n=l 



to the equations of motion and solving the resulting system of ordinary differential 
equations obtained from matching cos {ntot) terms: 



-' (r'^o)' 



<Po (00 - 1) (00 - 2) 

+ i (00 - 1) E (0n)' 
n 

+ i E 4>n4'p4'q iSn,±p±q) ■ 



(4.25) 



n,p,q 



r ^(r^^;)' 



n^w^ + 1 



+ 1 (00 - 1) E 0p0g (l^n.ipig) 

p,g 

+ 4 E 4>m.(t>n(f>q {Sn,±m±p±q) , 
m,p,q 



(4.26) 



where 6n,±p±q = {5n,+p+q + Sn,+p-q + K-p+q + K-p-q), and hkewise for dn,±m±p±q, 



but with nine terms. Equations 4.25 and 4.2f: can also be obtained by inserting 



ansatz 4.24 into the action and varying with respect to the 4>n ||5g] (see chapter g). 
This set of ODE's can be solved by "shooting", where the 0n(O) are the shooting 
parameters. Unfortunately, we were unable to construct a method that determined 



uj on its own; the best we could do was to solve equations 4.25 and 4.26 for a 
given CO that we measured from the dynamic solution. To easily compare the shoot- 
ing method to the dynamic data, we Fourier decomposed the dynamic data. This 
was done by taking the solution during the interval of time when the large period 



modulation disappears (1200 < t < 1800 for the oscillon in Fig. 4.9, for exam- 
ple) and constructing FFTs from the field at each gridpoint, rj. The amplitude of 
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each Fourier mode was obtained at each r, from the FFT using a window of 4096 
time-shces (bins). These amphtudes were then used to piece together the (/>„ (Fig. 



II2D 



Keeping only the first five modes, we compare the Fourier decomposed dy- 
namic data with the shooting solution (see Fig. 4.12| ). The value for to was de- 



termined from the dynamic solution and the shooting parameters (the </>„(0)) were 
varied so as to maximize rmaxj the distance before any mode diverged to (j) ^ ±00. 
The correspondence strongly suggests that the resonant solutions (ie. in the limit as 
ro -^ r^) observed in the dynamic simulations are indeed the periodic non-radiative 



oscillons obtained from ansatz ( 4.24 ) 



Looking at the three most dominant components of the power spectrum of 



(/)(0, t), Fig. 4.15 , we can see that during the period of no-modulation, the amplitude 
of each Fourier mode becomes constant. Although the plot is for the core amplitude, 
r = 0, this behavior holds for all r. This means that a "resonant" oscillon solution 
latches onto a non-radiative periodic solution as an intermediate attractor. Each 



resonance, however, does have a different critical solution; figure 4.11 shows distinct 
critical exponents for each resonance. This is reminiscent of the Type I critical 
phenomena studied in the massive Einstein-Klein-Gordon (EKG) model, |8[, Q, 
[PH], 1 18 1, where a family of (apparently) periodic solutions exist on the threshold 
of black hole formation. The behavior observed here is thus in contradistinction 



to Type I Einstein- Yang- Mills critical collapse, [15|, where the critical solution is 
found to be static and universal. In any case, instead of existing on the threshold of 
black hole formation, the critical oscillon solution is an oscillating time- dependent 
intermediate attractor on the threshold of having one more shape mode oscillation. 
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n = 2 



4 6 

r 



10 



Figure 4.12: Critical solution (l)n{r) (for n = 0,1,2,3,4) obtained from the Fourier de- 
composed dynarnic data (green x's) overlayed with (pnir) obtained by shooting equations 
4.25 and 4.26| (blue solid curves). The dynamic data was Fourier decomposed by forming 
a time series for each spatial gridpoint, r^, while the oscillon entered its non-radiative (no 
modulation) phase. FFT's were then constructed for each time series and the amplitude of 
each mode (the green x's) was measured. 
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Figure 4.13: Power spectra of the core amplitude, (/)(0,t), for the oscillons barely 
above and below the tq ~ 2.335 resonance. The power measured in each frequency 
regime slowly diminishes as the oscillon radiates away much of its energy until 
approximately t = 1100 where the oscillon enters a non-radiative state and all the 
components of the power spectrum become constant. 
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4.4.3 (Bounce) Windows to more Oscillons 

Lastly, we look at the region of the parameter space beyond rg ~ 5. The oscillons 
explored by Copeland, et al, were restricted to the domain of roughly 2 < tq < 5. It 
was even wondered if oscillons could form for larger initial bubble radii. We found 
that oscillons can form in this region and do so by a rather interesting mechanism. 
Again, going back to the 1+1 dimensional KK scattering of Campbell, et al, 
it is well known that the kink and antikink often "bounce" many times before either 
dispersing or falling into an (unstable) bound state. A bounce occurs when the 
kink and antikink reflect off one another, stop after a short distance, and recoUapse. 
We find that this happens in the (1+1) spherically symmetric case as well, where 
the unstable bound state is an oscillon. For larger rg, instead of the bubble wall 
remaining within r < 2.5 (as occurs for 2 < tq < 5), after reflecting off the origin 
the bubble wall can travel out to larger r (typically, 3 < r < 6), stop, and then 
recollapse, shedding away large amounts of energy in the process; we refer to these 
regions of parameter space as "bounce windows" . In the recollapse the simulation 
is effectively starting over with new initial data (albeit a different shape) but with 



much less energy (smaller tq). In Fig. 4.15 we see a lifetime profile and its resonances 
for a typical bounce window. 

4.5 The Klein-Gordon Equation in MIB Coordinates 
with ADWP {ID Critical Phenomena II) 

This chapter concludes with time evolution of spherically symmetric bubbles in the 



context of a slightly different action than the one introduced in section 4.1. The 
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Figure 4.14: Plot of (;/'(0,t) for rg = 7.25 displaying extremely nonlinear and 
unpredictable behavior during the "bouncing" phase (for t < 60), after which the 
field settles into a typical oscillon evolution. Once in the oscillon regime, the period 
is approximately uj ~ 4.5, and the first two modulations of the field can be seen 
(envelope maxima at i ~ 105 and t ~ 200). 
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Figure 4.15: Plot of oscillon lifetime versus initial radius of bubble for 4.22 < tq < 9. 
Although there seem to be no oscillons within 4.6 < tq < 6, it is clear that oscillons 
do exist for higher initial bubble radii. Furthermore, within these "windows" of 
parameter space that support oscillon formation, there also exist resonances similar 
to those observed throughout the 2 < rg < 4.6 region (see inset to observe resonances 
within 7.245 < ro < 7.257). 
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action used is the same as equation (^), 

S[cP]= Jd^x^\(^-^g^''V^cpV,<P-V{<P)^ , (4.27) 

except with an asymmetric potential, V{cf)) = VA{(t>) = -(f^ f 0^ — 4 (1 + 5) + 4j 

where again (j) = (j)(r,t), and (5 is a measure of the asymmetry of the potential. 

When (5 = 0, Va{4>) has the same shape as the SDWP potential of Chapter H, but is 

shifted so that the vacuum states are at i?i> = and (p = 2. With S ^ 0, the false and 

3 I 

true vacuum states are (f>F = and cfiT = - {I + 6) + - V 1 + 185 + 5^, respectively. 

The spherically symmetric MIB code is (trivially) modified by replacing the 
potential terms in the finite difference equations with the new potential, leaving b as 
a run-time parameter. Instead of performing extensive parameter space surveys with 



this potential like in section (4.4), a different type of phenomena is explored^. With 
the introduction of the asymmetric potential which has non-degenerate vacuum 
states, the threshold of expanding bubble formation can be examined. As discussed 
in section (p.l|), the initial bubble radius can be used to define a one-parameter 
family of initial data that transitions from bubble collapse (and eventual dispersal) to 
expanding bubble formation. As one familiar with critical phenomena might expect, 
there exists a time scaling law for the lifetime of the bubble lying on this threshold. 
The time scaling exponent is observed to be, 7 ~ 2.1, where T oc —7 In \a — a*\. 



^Brief parameter surveys were conducted and resonances were also observed, but no additional 
new resonance phenomena were observed. 

69 



0.5 



3 

< 

> 



■0.5 




Figure 4.16: Asymmetric double well potential Va((/>) 
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are at (ftp — and (/)t = — (1 
states, respectively). 



6) 



- V 1 + 18(5 + 6^ (the false vacuum and true vacuum 



70 




Figure 4.17: Plot of oscillon lifetime as a function of — ln|ro — rg| for the collapse (or 
expansion) of static bubble initial data using the ADWP. The lifetimes on the collapsing 
side of criticality are represented by green triangles, while the lifetimes on the expanding 
side of criticality are represented by red triangles. There is a single time scaling law with 
scaling exponent (slope) 7 « 2.1. 
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Chapter 5 



Axi-Symmetric Oscillon Dynamics 



This chapter discusses our investigations of the non-hnear Klein-Gordon model in 
axisynimetry. The equation of motion is written in 2-D MIB coordinates and a finite 
difference version of the equations is presented. The convergence properties of the 
code which solves the difference equations are shown, and evidence to support the 
validity of the MIB system as an effective absorbing outer-boundary is presented. 
A simple numerical technique is presented that allows for the arbitrary Lorentz 
boosting of spherically symmetric data. Finally, the code is used to collide two 
oscillons and a time-scaling phenomenon analogous to that seen in chapter 4.5 is 
observed. 



5.1 The Klein-Gordon Equation in Axi-Symmetric MIB 
Coordinates 

This section discusses the implementation of 2-D MIB coordinates in the solution of 
the nonlinear Klein-Gordon model with an asymmetric double-well (ADWP) poten- 
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tial. The model studied is the same ADWP model just used to study the threshold 
of expanding bubble formation in the previous chapter. The resultant PDE we must 
solve is still simply the Klein-Gordon equation 






but now with cp = (l){R, z, t). Defining 



U{R,z,t) = a6 (0 - /3^$R - r $, 

^R{R,Z,t) = dR(t>, 

^,{R,z,t) = 3,4), 



(5.1) 



(5.2) 
(5.3) 
(5.4) 



we have 



n 



2a- 



d 



d r:^ 



R /3^n + -<^n 
a 



R^ ,dV 

-^U-ab—j 

R dcf) 

ait \ab 

oz \ab J 



1 

ab 



oR 



n + /3«$R + /?^$, 



-|(^'"-^ 



(5.5) 
(5.6) 

(5.7) 
(5.8) 



for a, b, P^, P^, a, and R functions of R, z, and t, as defined in equations ( 3.27 ). As 
with the spherically symmetric case, these equations are similar in form to the Klein- 
Gordon equation coupled to gravity, except that, again, the "geometric variables" 
are known functions of our coordinates. 
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5.2 Finite Difference Equations 



Equations (|5.5D, ([5.6D, (p7?|), and ( |5.8| ) are solved using two- level second order (in 
both space and time) finite difference approximations on a static uniform mesh 
with N^ by N^ grid points in the R and z directions, respectively. The scale of 
discretization is set by AR, Az, and At = Amin(Ai?, Az)^. The field variables 
were first updated without dissipation everywhere (yielding the ' variables), and 
then dissipation was added were possible^ (yielding the ~ and final quantities). Using 
the difference operators from Table |3.2| , djiR = a, and dji = nR^^^djin, the interior 
{2 < i < Nfi — 1, 2 < j < Nz — 1) difference equations are: 



A,n», = 2rf" 



<,A,,(fl(„«n + i*„))+A,(rn + ^*.)]"^ 



n 

,ave ' ^' 



-fJ't 



R 



n-a6(/.((/>2-3(l + (5)0 + 2) I , (5.9) 



«j 



Ai(|.fi)J^. = ^--A^(^ln + /3«$^ + /3^^^)" , (5.10) 



*j 



1" .R. 



Atc&J.. = ^rA. H:n + /?''^fl + /3'^' , (5.11 



ab / jj 



1 \ " 

AtcPl^ = f,1--{ U + p'^^R + P'^') . (5.12) 

, 0,0 / j J 



For the inner R = boundary (i = 1, 2 < j < A^^ — 1), the conditions for regularity 
at the origin are applied for ^r and 11 (analogous to equations ( [4.15| ) and (4.16) for 
spherical symmetry): 



f^r 



A^n - ^n 



= (5.13) 

i,j 

^Although we usually take Az = AR 

^The sole purpose of this was to streamline the computer code. This allows the complete update 
to be execute using fewer loops. 
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H 



n+l 
id 



while evolving $2 and (j) a-s in the interior of the grid 



1 



^'1,3 \ab 



ab 



«j 



«j 




0. 



(5.14) 

(5.15) 
(5.16) 



For the outer boundaries, as with the spherically symmetric case, the equations 
used make very little difference since the physical position (i?,z) corresponding to 
the outermost gridpoint is moving out at nearly the speed of light. Therefore none 
of the outgoing field actually reaches the outer {R,z) boundary. Nevertheless, the 
typical massless scalar field outgoing boundary condition is imposed (on three edges 
and four corners): 
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for i = A^i?, 2 < j < iV, - 1, 
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1 
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^,J 



(5.17) 

(5.18) 

(5.19) 
(5.20) 

(5.21) 

(5.22) 

(5.23) 
(5.24) 



«j 
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for J = 1, 1 < i < Nr, 
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(5.25) 

(5.26) 

(5.27) 
(5.28) 



^,J 



for J = Rz, 1 < i < A^ij, 
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(5.29) 
(5.30) 

(5.31) 

(5.32) 



for i = 1, j = 1, 
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(5.34) 

(5.35) 
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*j 



for i = Nr, j = 1, 



^tli-lj 



f^t 



ra{^u + zAlu + n 

Vi?2 + z2 



«J 



(5.37) 
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(5.39) 
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for i = 1, j = iV^, and finally 
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(5.41) 

(5.42) 

(5.43) 
(5.44) 



«j 



for i = Nr, j = Nz- This condition is derived by transforming the spherically 
symmetric out-going boundary condition to axisymmetric coordinates. This is a 
reasonable approach to take since at large distances away from the collision, the 
emitted wave-fronts do becomes spherical. 

After each update using the evolution equations, dissipation is added inde- 
pendently in each spatial direction wherever the action of the appropriate dissipation 
operator is well-defined. Specifically, the fields are updated a second time (second 
updated quantities denoted by") using n'^^^ according to 



n: 



$ 



R 



(*•) 



n+1 
n+1 

n+1 



n 



n+1 



I dissfrn 



n+1 ,. / * \ n 



^■. 



n+1 



«J 



+ 4^"(C)" 



«J 



2 n+1 I ..diss 7 n 
^ij + ^^R 9i,j 



(5.45) 
(5.46) 
(5.47) 
(5.48) 
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where 3 < i < A^^ — 2 and 1 < J < -/Vz- Then the fields are updated a third time 
(third updated quantities have no accent) using /x^'*^*^ according to 

K^' = ni:;' + /xf^n^, (5.49) 

($ii)-;' = ($ij)"^^' + /if^(l'/j)"^. (5.50) 

($,)^+i = ($,)^^^. +/xf^(l>,)^^^. (5.51) 

'^"/' = 'Ai:;' + /^r "Ai:, (5.52) 

where 3 < j < N^ — 2 and 1 < i < A^^- The three update steps comprise one 
iteration. Typically, the whole process is repeated until the solutions at the (n + 1) 
time step converge to one part in lO^'^. 

5.3 Testing the 2-D code 

Although the 1-D MIB code worked very well, it was not obvious a priori that the 
same success would be achievable in a 2-D axisymmetric code. However, we find 
that the code is stable and passes all the desired tests; it is second-order convergent 
in general, it conserves energy to second-order in particular, and outgoing radiation 
is absorbed without refiection, so that there is no noticeable contamination of the 
interior solution (at least at the level of accuracy at which we work). 

Again, since equation ( |5.lD is a flatspace wave equation, obtaining a glob- 
ally conserved energy is straightforward. The spacetime admits a timelike Killing 
vector, f^, and therefore has a conserved current, J^ = f^T^u- A gaussian surface 
is constructed between two spacelike hypersurfaces within a cylindrical spatial do- 
main < R< i?cndcap, "^^endcap < z < ^endcap, with normals n^ = (±1,0,0,0) (the 
timelike normal to the hypersurfaces), n^ = (0, a^^, 0, 0) (the normal to the side of 
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the spatial gaussian cylinder), and n^ = (0,0, 0, ±6^^) (the normals to the top and 
bottom of the spatial gaussian cylinder). To obtain the conserved energy at a time, 
tf, the energy contained within the bubble, 

^endcap ^cndcap 



i?bubblc(t) = 2^ I J Rl^ + ^ + ^+V{ct>)]drdz, (5.53) 






U ^cndcap 

(where the integrand is evaluated at time t) is added to the total radiated energy, 
i?rad(t) = E^:^"^^^ + E^^"^^^ + i?Jd''''> where 

t ■^cndcap 

Ef:t'-{t) = 2^RJ J (^^yzdf (5.54) 

^ -^cndcap 

(where the integrand is evaluated at i? = i?endcap) and 

t -'T'cndcap 



^J^nd.ap (^) ^ 27r 11 R (^^) dR dt' (5.55) 



(where the integrand is evaluated at z = itzendcap)- The sum, -Btotai = -E'bubbic+-E'rad) 
remains conserved to within a percent^ at 201 x 401 gridpoints (see figure 5.1 for 



energy as a function of time, see figure ^^ to see more clearly the energy being 
conserved as per second order convergence). In addition to conserving energy (to 
second order), the code displays second order convergence in all the evolved field 
variables (see figure |5^ ) . However, convergent energy conservation and convergent 
field variables do not imply that the MIB system is absorbing the outgoing radiation 
properly. 

To verify that there is no significant refiected radiation off the outer boundary 
that contaminates the solution (exactly following the methods described in section 
|4.3| ) the MIB solution is compared to an ideal large-grid solution by taking the 

^ Using i?niax = 25 and 2max/min = ±25. 
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Time 



Figure 5.1: Energy conservation of 2D code at four discretizations, 41 x 81, 81 x 161, 
161 X 321, and 321 x 641 gridpoints (8h, 4h, 2h, and h, respectively). The code conserves 
energy approximately to one part in 13, 50, 200, and 800, for levels 8h, 4h, 2h, and h, 
respectively, thus displaying energy conservation consistent with a second-order convergent 
code. 
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Figure 5.2: Convergence factors, {E4h — E2h) / {E2h — Eh) in red and 

{E^h — Eih) / {Eih — E2h) in g reen , for the conserved energy. As with the conver- 
gence tests in chapter Q (figure O, in particular) a convergence factor of four indicates 
second-order convergence. The code clearly conserves energy to second order. 
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Figure 5.3: Convergence of 2D code for the field at four discretizations, 41 x 81, 81 x 161, 
161 X 321, and 321 x 641 gridpoints (8/i, Ah, 2h, and h, respectively). The convergence factors 
for levels 8h, Ah, and 2h and levels Ah, 2h, and h are shown in red and green, respectively. 
The value of roughly four observed for each convergence factor indicates that the code is 
indeed second-order. 
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f2-norm of the difference at every point, 

1/2 



El«uf I (5-56) 



MN ^ 



where Ojj are the components of an M x A^ matrix, A. This norm is compared to 
the ^2-iiorm of the difference between the large-grid solution and a solution using an 
out-going boundary condition known not to work well (again, the massless outgoing 



boundary condition, OBC). Figure 5.4 shows that after two crossing times, the 
OBC solution becomes contaminated and the solution error increases dramatically 
(approximately three orders of magnitude!) while the MIB solution always remains 
around or below 10"^ and shows no dramatic increase in solution error indicative 
of contamination. 

Just like with the 1-D code, since the 2-D MIB code conserves energy quadrat- 
ically, has quadratically convergent field variables, and is in agreement with the 



large-grid solution; it is an acceptable means of solving equations (^^),(^^,(5.7) 



and ( |5.8| ), while being dramatically more computationally efficient than either dy- 
namical grid methods or large-grid methods. The computational demand for the 
MIB system grows linearly with the oscillon lifetime, while for dynamical or large- 
grid methods, the computational demand grows as the cube of the oscillon lifetime^. 

5.4 Boosting the Spherically Symmetric Oscillons as 
Initial Data 

The main (certainly most fun) reason for creating an axisymmetric code to solve 
the KG equation was to investigate what happens when two oscillons collide. 



Assuming that the outgoing radiation demands the grid to grow in both the R and z direction. 
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Figure 5.4: Plot comparing the OBC (dashed red) and MIB (solid green) solutions to 
an "ideal" solution. The OBC solution is obtained from using a massless outgoing bound- 
ary condition, the MIB solution is obtained by solving the system in axisymmetric MIB 
coordinates, and the ideal solution is obtained by evolving the solution in standard {R,z,t) 
coordinates on a grid large enough to ensure no reflection off the outer boundaries. The 
error estimates are obtained from the ^2-norm of the difference between the trial solutions 
(OBC or MIB) and the ideal solution, ||0 — 0idcai||F. Contamination of the OBC solution 
is observed at two crossing times, t « 40, where the error estimate increases almost three 
orders of magnitude. 
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The first step in most dynamical calculations is to generate meaningful initial 
data. Since there is no "closed-form" oscillon solution, it is actually non-trivial to 
generate "boosted" oscillon data. There are many approximate initial data config- 
urations that give rise to a gaussian-shaped bubble that can move at the desired 
velocity. For example, "boosting" 

(/)(i?, z, i) = 00 exp 2 2 H^"^ ("^^ (5.57) 

\ '^R '^z ) 

by applying the coordinate transformation for a Lorentz boost along the z-axis (using 

( R"^ (-i(z -vt))'^\ 
(t>{R, z, t) = 00 exp ^ - ^-^^ — 5-^^ sin {co-/ {t - vz)) , (5.58) 

V ^fl '^i J 

usually results in an oscillon that moves at roughly the desired velocity^. The 
problem, however, is that even with the best possible choices of the free parameters 
{(po, uj, aji, Gz) the approximate oscillon ( |5.57 ) is not a solution to the (Lorentz 



invariant) Klein-Gordon equation. Therefore, trying to generate initial data in this 
manner often results in dramatically different global behavior depending on the 
"boost velocity" , v (ie. behavior that is not Lorentz invariant and can be as dramatic 
as leaving the universe in different vacuum states!). The goal is to be able to compare 
oscillons being collided at different velocities while knowing exactly how each oscillon 
behaves in its own rest-frame. 

We choose to numerically evolve the oscillon in its rest frame (which can be 
done efficiently in spherical symmetry) and then numerically boost the solution. The 
single boosted oscillon then, of course, retains all of its Lorentz invariant properties 



^Vacuum assumed to be at </> = 0. 
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regardless of boost velocity. The collision initial data is constructed from the super- 
position of two boosted oscillons (one at some z < boosted in the z~^ direction, 
and the other at some z > boosted in the z^ direction). 

The equations for the transformation between the boosted (tilde) and rest 
(non-tilde) frames are just the equations for a Lorentz boost along the z-axis: 

t = J (t — vz) t = jit + vz) 

z = ^ (z — vt) z = ^ (z + vi) 

(5.59) 

R ^ R R ^ R 



where again, 7 = :. Since is a scalar field, it transforms trivially as 

V 1 — v'^ 

4>ii, R, 6i, z) = ((7 (t - vz) ,R,e,j{z- vt)) , (5.60) 



while the other field variables (being derivatives of the scalar field) transform as 



forms d^'(p = ^ ^^, d^(j) which gives 



n = 711 -I- 'yv^2 



jU+ J^ ^r. (5.61) 

VR^ + z^ 



where $,■ = dr(j) and is obtained from the spherically symmetric evolution. (Re- 



member that r = a/x^ + y^ + z^ and R = J x"^ + y"^.) 

Equations ( 5.59| ) imply that to obtain data for t = 0, Rm\n < -R < -Rmax, 



and Zmin ^ z < imaxj the rest frame solution must be known over the spacetime 
domain ^ivz^^i < t < 7UZmax, ^min < R < -Rmax, and 7Zmin < z < 75max- Time 
symmetric initial data was used and the z domain was chosen such that Zmin = 
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Figure 5.5: A schematic (for a sample constant-i? slice) displaying how the {R,z,io) 
boosted data is obtained by interpolating from rest-frame {R, z,t) dynamic data. Data 
points that do not match up exactly between the tilde and non-tilde frames are linearly 
interpolated in both space and time (note that figure is only a schematic and not drawn to 
scale, Az and At were always taken to be very small compared to Az). The blue squares 
are gridpoints lying on t-constant hypersurfaces in the rest-frame. Red circles are gridpoints 
lying on the t = boosted- frame hypersurface (the initial data being obtained) . To obtain 
boosted data for Zmin "£ z < Zmax at t = 0, the rest-frame radius must be known over the 
domain jvzmin < ^ < jvzmax and jZmin < z < J Zmax- Since R and R are orthogonal to 
the boost, the field remains unaffected in that direction. 



which allowed the time domain for the spherical evolution to be restricted to 



< t < jv Zmax (see figure x5 for a schematic relating a "new" t = slice to the 
"old" spacetime {i,z) domain, for a typical R = constant slice). 

A boosted oscillon can then be obtained by time evolving the desired (rest 
frame) oscillon throughout the necessary spacetime domain and transforming the 
field variables using equations ( 5.60| , 5.61 ). Since very few of the (t, R, z) gridpoints 
match exactly any gridpoint from the spherically symmetric (i,r) grid, linear inter- 
polation in both space and time was used (the 1-D grid that was interpolated from 
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is at a much finer resolution than the 2-D grid). A schematic of an constant-i? slice 
can be seen in figure ^^. Finally, collision initial data is obtained by forming the 
superposition of two boosted oscillons (obviously, boosted at one another.). 

5.4.1 Testing the Numerically Boosted Initial Data 

The "numerical boosting" of the spherically symmetric oscillon data was tested by 
monitoring the Lorentz invariant behavior of expanding bubble formation. For a 
spherically symmetric gaussian pulse, instantaneously at rest, 

(j){r,0) = cj)F + {'pT-'pF)expl-^\ (5.62) 

(/)(r,0) = (5.63) 

so </)(c«,t) = (pF and (l){0,t) = (pT, there will be a critical "bubble radius", a*, for 
which time evolution of ar > cr* will lead to an expanding bubble that converts 
false vacuum to true vacuum everywhere. For ar < cr* the time evolution will 
eventually result in dispersal of the field (possibly after an oscillon stage) . Since a* 
is defined to be the rest-frame critical radius for expanding bubble formation, this 
number should not change and the boosting of the critical field configuration should 
always lie on the threshold of expanding bubble formation. Figure ( |5.6D shows the 
critical bubble radius as a function of boost velocity for three resolutions (81 x 161, 
161 X 321, 321 X 641); the critical bubble radius remains constant to within one- 
quarter percent at 161 x 321 for v < 0.8, and to within one percent at v = 0.85. The 
three resolutions show convergent properties (ie. the graph of a* vs. boost velocity 
approaches a constant as /i ^ 0) and suggest that the deviation from the expected 
value of a* at high boost velocity is due to the steep gradients which primarily result 
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Figure 5.6: Plot of the Critical Initial Bubble Radius versus Boost Velocity for three 
discretizations (81 x 161 in red triangle, 161 x 321 in green squares, and 321 x 641 in 
blue circles). Since the initial bubble radius parameter in the code is the rest-frame 
initial radius, the Lorentz invariant behavior of contraction or expansion (which 
governs the threshold value of tq) should not change as a function of boost velocity. 
The "ideal" result is the (dashed cyan) horizontal line at rg ~ 3.34 (determined 
spherically symmetric collapse at N^ = 2049). The deviation at large boost velocity 
is due to the increasing effects of length contraction of the bubble. 
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from the length contraction of the bubble. 

5.5 Collision and Endstate Detection 

In the spherically symmetric ADWP collapse, there were two possible endstates: 
expanding bubble formation or bubble collapse (usually resulting in the formation of 
an oscillon). However, in axisymmetric collision simulations, there are four possible 
endstates that we consider^: expanding bubble formation, annihilation, soliton-like 
transmission, and coalescence. 

Expanding bubble formation occurs when the whole spacetime is converted 
from the (unstable) false vacuum to the (stable) true vacuum. This can occur 
trivially when the bubbles that are being boosted at one another each have initial 
rest-frame radii above the critical radius for expanding bubble formation (for 6 = 
0.1, a* ^ 3.3); this is the case where each bubble would lead to an expanding 
bubble independently. However, a more interesting case involves the collision of two 
oscillons, each with a rest-frame radius less than the critical radius for expanding 
bubble formation. In this case, each oscillon independently would not lead to an 
expanding bubble, but together they overcome the potential barrier of Va_{4>) and a 



false-to-true phase transition ensues (see figure 5^ for heuristic explanation). The 
code detects bubble formation by computing the area {J dR dz, since we work in 
two spatial dimensions) that is in the true-vacuum and stopping the evolution after 
some empirically determined threshold is reached. In practice, it is quite easy to 

see where the "surface tension" in the bubble wall can no longer compete with the 

®Of course, there could be more ways of classifying them into more endstates, but we choose 
four. 
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Figure 5.7: Heuristic representation of what happens energetically when two oscillons col- 
lide to form an expanding true vacuum bubble. Individually, the two bubble configurations 
(represented by the two balls oscillating about the false vacuum) are far from one another 
and do not individually have enough energy to overcome the potential barrier; if left alone, 
they would eventually disperse, leaving the spacetime in the false vacuum. However, when 
the two oscillons collide, their combined energy is enough overcome the potential barrier 
and form an expanding bubble that converts the spacetime to true vacuum. A time series 
of an actual collision of this type can be seen in figure 5.8. 
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field's "attraction" to the true-vacuumQ; the detection threshold is then taken to be 
safely above this point. A time series of snapshots of an actual 2-D collision that 
results in the formation of an expanding bubble can be seen in figure |5] 

The other three endstates being considered do not form expanding bubbles 
and therefore leave the spacetime in the false vacuum. Two of these three endstates, 
annihilation and soliton-like transmission are forms of dispersal, while coalescing 
oscillons remain within the computational domain. To help distinguish between 
these endstates, three second moments are considered^: 

/ R^p dR dz 

Ml{t) = J— (5.64) 

/ p dR dz 

z p dR dz 

Mlit) = J— (5.65) 

p dR dz 

R^ + z^) p dR dz 



where 



M;(t) = ^^— — (5.66) 

/ p dR dz 



1 / n2 $2 ^2\ 

p is the energy density due to the time derivative and gradients of the field, neglecting 
the potential terms. Including the potential term would lead to large contributions 
to the moments from outside the bubble, whereas the goal here is to know the the 
location of the surface of the bubbles. For all of the collisions discussed here, the 
initial data are symmetric about the z = plane, and hence the collision also occurs 
at z = 0. 



^We are also guided by critical radius measured from the 1-D ADWP results. 
*Note that these are clearly not physical moments, but they work well to distinguish between 
the various endstates of the oscillon-oscillon collisions. 
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Vb 


(^0 


A\v\ 


AT 


0.6 


2.15 


+0.02 


-1.7 


0.6 


2.675 


-0.41 


-15.8 


0.85 


2.00 


-0.02 


-0.6 


0.85 


3.325 


-0.17 


-0.5 



Table 5.1: Table of change in the osciUon velocity (magnitude), A|w|, and time "lag", AT, 
for four sample soliton-like collisions. Oscillons (and solitons) that interact nonlinearly are 
usually characterized by a time shift or change in velocity arising from coupling between 
modes. If the oscillons interacted linearly, A|w| and AT would be zero. The oscillons 
explored here almost always slow down or stay at roughly the same speed and the time lag 
(measured by extrapolating back to z = 0) almost always is negative. This is contrary to 
what is usually seen in (1+1) dimensional soliton interactions where AT > 0. 



Annihilation occurs when two oscillons collide and interact in such a way that 
the field is no longer localized and all the radiation disperses. The code determines 
this to be the endstate when both M^ and M^ rise above an empirically chosen 
threshold. This was typically taken to be around 100 which implies most of the 
matter is around r = yB? + z^ = 10), which, in turn, is indicative of dispersal for 
oscillons with typical radii around r ~ 3. A time series of a collision resulting in 
annihilation can be seen in figure p.9| . 

Soliton-like Transmission occurs when two oscillons collide, interact, and 
then pass through (or reflect off) one another while each oscillon stays localized. 
Again, this is easily determined by monitoring the second moments. Since the os- 
cillons pass through one another, the z moment will get large as they leave the 
computational domain. However, since each oscillon remains localized, the R mo- 
ment will stay small (on the order of the typical oscillon radius of r ~ 3). A time 



series of a collision resulting in soliton-like transmission can be seen in figure 5.11. 

Lastly, Coalescence occurs when the two oscillons collide, interact, radiate 
away some of their energy, and survive to form a single oscillon. Coalescence is 
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(Area)rpj,^g vacuum 


Ml 


M'i 


Expanding Bubble 


large 


- 


- 


Annihilation 


small 


large 


large 


Soliton-like 


small 


small 


large 


Coalescence 


small 


small 


small 



Table 5.2: Logic table for endstate classification. (Area)rpj.^^ vacuum ^^ ^^^ "area" (/ dz dR) 
of space that is in the true vacuum. M|, and M^ are the second R and z moments, respec- 
tively, of the energy density-like function, p. The expanding bubbles are easily separated 
from the other endstates; once the area of space at the true vacuum reaches a given thresh- 
old, the collision is flagged as forming an expanding bubble. Annihilation and soliton-like 
transmission both occur when the field disperses (M^ becomes large) and are differentiated 
by the mass moment in the radial direction; if M^ is small, the field is localized indicative 
of transmission, whereas if M^ is large, the field has clearly dispersed. 



the default endstate for the code's detection algorithm since it is what happens 
if an expanding bubble does not form or the field does not disperse in the form 
of annihilation or soliton-like transmission. Unfortunately, simulations with this 
endstate are the most computationally intensive because it is hard to know how 
long to wait before deciding that the solution will not disperse or form an expanding 
bubble (particularly on the threshold of expanding bubble formation) . A time series 



of an oscillon-oscillon collision that results in coalescence can be seen in figure 5.12. 



The endstate classification logic is summed up in Table 5.2. 



5.6 Parameter Space Survey 



In the spherically symmetric ADWP evolutions, a natural parameter to adjust in 
the exploration of the model was the initial width of the collapsing bubble, ar of 
equation ( |5.62 ). This one dimensional parameter space is extended here by adding 
the boost velocity, Ufe, as a second parameter. Although the initial data is still 
highly symmetric (the colliding oscillons are the same size and oscillate in phase), 
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there is still interesting structure to the parameter space. In figure 5.13 we see that 
for a large area of parameter space, there are collisions where bubbles that would 
otherwise have (eventually) dispersed actually combine to form an expanding true 
vacuum bubble. 

5.7 Threshold of Expanding Bubble Formation in 2D 
Collisions 

{2D Critical Phenomena) 

At the end of chapter ^, a time scaling law was presented for solutions on the 
threshold of expanding bubble formation for spherically symmetric bubble collapse. 
With the 2D parameter space surveys presented in section |5.6| (for collisions) , there 
also is a boundary separating expanding bubble formation from solutions that leave 
the space in the false vacuum. The natural question is then: Does a similar scaling 
law exist, and if so, does it have the same critical exponent? 

The threshold of expanding bubble formation was explored first by fixing 
the boost velocity, v^, and then bisecting between two values of the initial rest 
frame radius, ar (one that forms an expanding vacuum bubble and one that does 
not). Next, the threshold was explored along a different direction in parameter 
space by fixing ar and bisecting between two values of Vb. Time scaling relations 
were observed for each bisection and the critical exponents were observed to be 
quite similar (7 ~ 2.2, see table ^.3]) . The locations in parameter space of these 



critical points can be seen by the crosses in figure 5.13. The similarity between 



the critical exponents for the axisymmetric collision of two gaussian bubbles and 
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< 


AfT 


vl 


Av*, 


7 


2.620314246044920 


~ 10-^4 


0.4 


fixed 


2.18 


2.824162086503015 


~ 10-14 


0.5 


fixed 


2.15 


2.813209312722570 


~ 10-14 


0.6 


fixed 


2.19 


3.293020631075245 


~ 10-14 


0.6 


fixed 


2.15 


3.281889627766725 


~ 10-14 


0.7 


fixed 


2.15 


3.270433308449535 


~ 10-14 


0.75 


fixed 


2.15 


3.278874383866785 


~ 10-14 


0.8 


fixed 


2.17 


2.5 


fixed 


0.292882455348505 


~ 10-14 


2.18 


2.6 


fixed 


0.206552206744452 


~ 10-14 


2.18 


2.65 


fixed 


0.415868755640645 


~ 10-14 


2.17 


2.7 


fixed 


0.337408336301504 


~ 10-14 


2.17 


3.0 


fixed 


0.658702864946936 


~ 10-14 


2.20 


3.0 


fixed 


0.327480589060684 


~ 10-14 


2.17 


3.2 


fixed 


0.546131322091751 


~ 10-14 


2.19 



Table 5.3: Table of critical exponents and threshold parameter values for a two- 
dimensional parameter space survey of axisymmetric collisions. The entire survey can be 
seen in figure p. 12 , where the black x's denote the threshold values in this table. The points 
were chosen to span the space as much as possible and the bisection search was performed in 
both directions of the parameter space; nevertheless the critical exponents are remarkably 



similar, suggesting a universal scaling law (see figure 5.15) 



the spherically symmetric collapse of a single gaussian bubble suggests that the 
dominant unstable mode observed with the axisymmetric critical solution might be 
spherically symmetric. 
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(a) 



(b) 





(c) 



(d) 



Figure 5.8: A times series showing the collision of two oscillons forming an expanding 
true vacuum bubble. The initial frame (a) shows the two oscillons boosted toward one 
another; the fields interpolate between the vacua so that their peaks are at the true vacuum 
while away from the oscillons is the false vacuum. Frames (b) and (c) show the nonlinear 
interactions during the collision. Frame (d) shows the beginning of an expanding bubble; 
the steep field gradient is the bubble wall. Inside the bubble wall is the true vacuum and 
the bubble wall will continue to move outward (eventually reaching the speed of light). 
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(a) 



(b) 





(c) 



(d) 



Figure 5.9: A times series showing the collision of two oscillons that annihilate each 
other. The initial frame (a) shows the two oscillons boosted toward one another; the fields 
interpolate between the vacua so that their peaks are at the true vacuum while away from the 
oscillons is the false vacuum. Frames (b) and (c) show the nonlinear interactions during the 
collision. Frame (d) shows the endstate consisting of almost entirely of outgoing radiation. 
The spacetime is left in the false vacuum. 
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Time 




(a) Vb = 0.6, (To = 2.15 



(b) Vb = 0.85, CTo = 2.0 





Time 
(c) Vb = 0.6, ao = 2.675 



Time 
(d) Vb = 0.85, ao = 3.325 



Figure 5.10: Plots of the location, (Af^)^/^, of the oscillon located in the "upper" (z > 0) 
half of the computational domain as a function of time. In either plot, the green line is the 
best fit for the ingoing position as a function of time, the blue line is the best fit for the 
outgoing position, and the cyan line is the path that the oscillon would have taken had it 
not interacted (nonlinearly) with the other oscillon (direct reflection across z = 0). During 
the transmission the velocities (slopes) can increase or decrease while the outgoing pat h was 
always seen to be shifted backward in time (x- intercept always moved left). Table ^.l| shows 
the measured time lags and changes in velocity for these four evolutions. 
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(a) 



(b) 





(c) 



(d) 



Figure 5.11: A times series showing the cohision of two oscillons that pass through (or 
reflect off) one another. The initial frame (a) shows the two oscillons boosted toward one 
another; the fields interpolate between the vacua so that their peaks are at the true vacuum 
while away from the oscillons is the false vacuum. Frames (b) and (c) show the nonlinear 
interactions during the collision. Frame (d) shows the endstate consisting of two oscillons 
moving outward along the z-axis. The spacetime is left in the false vacuum. 
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(a) 



(b) 





(c) (d) 

Figure 5.12: A times series showing the collision of two oscillons that coalesce. The initial 
frame (a) shows the two oscillons boosted toward one another; the fields interpolate between 
the vacua so that their peaks are at the true vacuum while away from the oscillons is the 
false vacuum. Frames (b) and (c) show the nonlinear interactions during the collision, where 
roughly half the total energy of the system is converted into radiation. Frame (d) shows the 
endstate consisting of one remaining oscillon, oscillating on top of the false vacuum. 
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Figure 5.13: Two dimensional, CTQ-Vb, parameter space survey of axisymmetric 
bubble collisions (colored according to the endstate of the system); the vertical 
axis measures the initial rest-frame radius of the bubble, ctq, while the horizontal 
axis measures the velocity that each oscillon is boosted at, Vb- The red region 
corresponds to expanding vacuum bubble formation, the cyan region to coalescence, 
the blue region to soliton-like transmission, and the green to annihilation. The 
black curve is the line representing the threshold for expanding bubble formation for 
individual boosted bubbles (bubbles above this line would form expanding vacuum 
bubbles even if not colliding). The blacl^O?'s correspond to points on the threshold 
of expanding bubble formation for which critical phenomena was explored (see table 
p.3|). Each collision was performed on a 161x321 grid. 
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(a) 81 X 161, h = 0.25 



(b) 114 X 227, h « 0.177 





Boost Velocity 
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(c) 161 X 321, h = 0.125 



(d) 227 X 453, h « 0.088 



Figure 5.14: (Jq-vi, parameter space surveys at four different discretizations (Ar = Az ~ 
h). Although no particular convergence factor was measured, the shapes of the colored 
regions appear to be converging to those of figure (d) , at least well enough to support the 
existence of both the large oval-shaped redlft^ion centered around vb ~ 0.5, (Tq ~ 3 and 
the side- lobes that are present for smaller boost velocities (to the left of the main oval). 
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Figure 5.15: Sample Time versus — In \ao — Cq I plot for tlie tlireshold of expanding bubble 
formation. The points plotted are for gq > CTq and therefore all produce expanding bubbles. 
The scaling exponent (slope) was measured to be 7 « 2.15. All of the th resh old parameters 
explored produced similar scaling exponents and are tabulated in figure 5.3. 
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Chapter 6 



Trial Function Approach to 
Critical Solutions 



This chapter discusses both exact and approximate critical (non-radiative) oscillon 
solutions obtained through a trial function approach. Differential equations ( 4.25| ) 



and ( [4.26| ) are rederived using a very general trial function ansatz. By assuming 
more constrained gaussian and hyperbolic secant trial functions, approximate critical 
oscillon solutions are also obtained. 

6.1 Trial Functions &; Variational Approach 



In section 4.4.2 , the ordinary differential equations ( 4.25 ) and ( 4.26 ) were derived 



by inserting the ansatz 

oo 

0(r, t)=-l + Mr) + J2 Mr) cos (nujt) (6.1) 



n=l 
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into the equation of motion ( [4 .21 ) and matching cos{nujt) terms. The same set of 
differential equations can be obtained by inserting ( |6.lD into the action 

S = I (fx C 



(c o\ 



and integrating over time to obtain a dimensionally reduced action^. The reduced 
action is then varied with respect to the individual modes, (/>o and 0„, to obtain the 
(ordinary differential) equations of motion. 



Ansatz ( 4.24J6.1 ) is actually obtained by starting with the most general pos- 
sible ansatz 

oo 
(/>(r, t) = -l + (t)o{r) + Y, ^n{r) sin {nujt + On) (6.3) 

n=l 

and imposing boundary conditions that match the phases of the each mode and that 



demand that the maximum amplitude of each mode is equal to An (0) . The reduced 



action obtained by integrating ( |6.2| ) and (6J) over time is 

/•oo rT 

S = An dr C{r,t) 

Jr=0 Jt=0 

" /-oo f 1 °° o 1 °° 

2f--{6cl>l-12cPo + 4)Y^ 



UJ Jr=0 



71=1 

- (00 - 1) Yl (t'n'Pp^q iK±p±q) 
n,p,q=l 
-. oo 

-T7 z2 4>n4>m(t>p4>q {^n,±m.±p±q) 
n,m,p,q=l 



4^i /-oo 

/ dr C{r) 

U) Jr=0 



^Of course, the second line in ( |6.2| ) is also dimensionally reduced by integrating over dfl, which 
is trivial in spherical symmetry. The more interesting case arises when imposing the (discrete) 
periodic boundary condition in time. 
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where again^ the pluses and minuses sum over every permutation (ie. Sn,±p±q = 
{6n,+p+q + Sn,+p-q + Sn,-p+g + Sn,-p-q))- Unfortunately, ansatz ( |6.1[) is so general 
that the spatial integration cannot be done in closed-form. However, one can still 
vary the (partially) reduced action with respect to (po and (/>„ to obtain the usual 
Euler-Lagrange equations, 

f dCjr) \ dC{r) 

^ ( dC{r) \ dCir) , , 

to obtain a set of differential equations which (again) are 

1 Ir^'iAl 






r^<t>'o)' = </.o (<^o - 1) (<^o - 2) 



+i ('Ao - 1) E (</'«)' (6.7) 

n 

+ \ E 4>n(kp(t>q {^n,±p±q) , 
n,p,q 



{r^<P'J = (3icP,-lf-{nW + l] 



1 f^2j,/ 



+ 1 {(t>0 - 1) E (t>p<Pq {K,±p±q) (6.8) 



+ 4 E (t>m<t>n(t>q {K,±m±p±q) , 
m,p,q 



respectively. The solution to equations ( |6.7| ) and ( |6.8|) for uj = 1.38 can be seen 
in Figure 4.12| . The solution matches extremely well to the Fourier decomposed 



dynamic data. Again, this is not surprising because the only assumption made 
restricting the ansatz that could prevent solutions to ( |6.7D and ( |6.^ ) from being 
solutions to the true PDF equations of motion was that the solution is periodic, a 
property which is observed empirically for the resonant oscillons. 



/" /7r\ 

^ We also made judicious use of / cos{mujt) cos{nLot) — { — ] Sm n, for m,n > 1. This is why 

Jt=o ^"^^ 

the ansatz was split up into the (jto and (j}„ modes, since the above integral is 2n/uj for m = n = 0, 

not Tv/ijSm.n as for m, n > 1. 
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6.2 Approximate Oscillons From Constrained Ansatz 

Although the resonant non-radiative solutions found in section ^T^ match the dy- 
namic data extremely well, they are still obtained from solving ordinary differential 
equations. With a well chosen and more constrained trial function, both the time 
and space integrations in the action can be evaluated in closed-form. Variational 
methods are then applied to obtain a set of algebraic equations that can be solved 
to obtain the best possible values for the free parameters of the trial function. The 
ansatz is taken to be 

(j)(r,t) = -l + A f(-jsm{Lot + d), (6.9) 

where a is related to the width of the field configuration and A is the amplitude of 
the oscillations. Equation ( |6.9D is subjected to the following boundary conditions: 

0(r,O) = -1 (6.10) 

Hr,7r/{2u;)) = -1 + Af (^^^ (6.11) 

which results in a choice of phase, 9 = 0, and sets A equal to the maximum amplitude 
of the oscillating mode. This yields 

(t){r,t) = -l + A f(-jsm{ujt) (6.12) 



for the general form of the trial function. Equation ( 6.12| ) is inserted into the action 
and can be written in terms of definite integrals. 



' - '' f'^ir' "'■'{¥' -I*" -\i*'-'y 



IT [2A^a^uJC2gi - uJ "^ {2A^acig2 + 4.A^a^C2g2 - 4^^0-^0353 + A^a^c^g^j | 

(6.13) 
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where the q are dimensionless constants 

/•OO /"OO 

ci = / du u {duf{u)) C2 = / du u f{u) 

Jo Jo (g_^4) 

C3 = / du u f{u) C4 = / du u f{u) 

Jo Jo 

obtained by taking u = r/a and the gi are dimensionless constants 

91 = dr cos (r) = - 

JO 4 



rV2 _ , vr 

4 



r' 2 

92 = dr sin (r) 

•^° , " (6.15) 

93 = dr sur^T) = — 

Jo o 

rTT/2 3yp 

54 = / dr sin'^(r) = — 

Jo 16 

obtained by taking t = ut. The reduced action ( |6.13| ) will be the same even if 
sin(u;t) in equation (6.12) is replaced with an arbitrary function, h{ojt). However, 
the constants, 51, §2, 93, and 54, will have to be recomputed according to 



(6.16) 



91 = r dr {drh{T)f 52 = r dr h{Tf 

Jto Jtq 

93 ^ rdrhirf 54 = rdrhirf 

Jto Jto 

and the limits of integration need to be chosen appropriately. 

Now that the action is completely reduced, varying with respect to the free 
parameters yields a system of algebraic equations, 

(6.17) 
(6.18) 

that can be solved for A and a^, 



dS 
da 


= 


as 
dA 


= 0, 



6C353 ± JSGcjgl - I6C2C494 (2^2 - 51^^) 
A± = "^ (6.19) 

20454 



2 -2ci52 



""^ 3 (C454^^ - 4c353^± + 2c2 (252 - 51^')) ■ ^^'^^^ 
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The reader may be wondering why, since the action is also a function of w, 
that the equation — — = was not added to the system of equations, ( |6.17 ) and 



( |6.18| ). Unfortunately, the trial function method does not work well to determine 
io. Looking at typical plots of the action^ as a function of the amplitude and width 



(figures 6.1 and |6.2| , respectively), it is apparent that there are no (non-trivial) 



points where — — = 0. Forcing this equation to be solved in conjunction with (|6.17| ) 
and ( |6.18| ) results in unphysical complex solutions. In the following subsections, we 
choose three forms for f{r/a) and look at the amplitude and widths as a function 
of a;. 

6.2.1 Gaussian Trial Functions 

In this subsection the field configuration is assumed to have a gaussian shape (as 
assumed in [|2^), 



Now that f{r/a) has been chosen, the integrals ( 6.16| ) can now be performed: 



3 /vF 1 [¥ 



3 V 2 8 V 2 (g_22) 

1 f¥ 1 



^^ = I2V3 '' = 32^"- 



Using equations (|6.19| ) and ( 6.2C| ), A± and a"^ can be plotted as a function of 



CO (figure S.3). Unfortunately, for cu < 0.437 the discriminant in equation ( |6.19 ) is 
negative and gives complex solutions. Furthermore, where A- is real cr^ is negative, 
resulting in an imaginary oscillon width! However, there is a part of the "+" root 

(0.62 < a; < 1.41) that has positive o"^ and can therefore support physical solutions. 

f{r/ij) is taken to be a gaussian for these figures, but the property holds for the sech and sech 
functions also discussed below. 
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omega 



Figure 6.1: Plot of the reduced action as a function of the oscillon amphtude, A, and the 
frequency, u; (while setting a — 3.6). Looking along a;-constant slices, it is clear that there 

dS 
are many points that satisfy -— - — 0, whereas when looking at A-constant slices, there are 

dS 
no points that satisfy — — = (except for the trivial A — solution). 

OLU 
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Figure 6.2: Plot of the reduced action as a function of the osciUon width, ct, and the 
frequency, w (while setting A = 0.62). Looking along w-constant slices, it is clear that there 

dS 
are many points that satisfy — — = 0, whereas when looking at cr-constant slices, there are 

da 
dS 
no points that satisfy — — = (except for the trivial a = solution). 
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Parameter 


Dynamic 
Solution 


A+ = Ad 


a+ = (TD 


A 
a 

UJ 


0.62 
3.6 

1.4 


0.62 
2.9 
1.26 


0.39 
3.6 

1.32 



Table 6.1: Table of free parameter values {A, a, and lo) determined using a gaussian trial 
function (third and fourth columns) compared to the values estimated from the spherically 
symmetric dynamics (second column); the D subscript indicates the value was obtained from 
the dynamic solution. The third column contains results obtained using the trial function 
approach where A^ equals the amplitude of the n = 1 mode dynamic solution; the fourth 
column contains results where cr+ equals the width of the n = 1 mode dynamic solution. 
The data from the third and fourth columns occur around the x's in figure |6.3|. 



Table 6.1 compares sample solutions (located near the x's in figure 6.3) to estimates 
obtained from the dynamic solution. 



6.2.2 Hyperbolic Secant Trial Functions 



In this subsection the field configuration is assumed to have a hyperbolic secant 
shape, 

/(^^)=sech(^^). (6.23) 

With the new choice of f{r/a), the integrals ( |6.16| ) are now reevaluated: 

1 2 1 1 2 

Cl = IT ^ Co = TT 

36 3 ' 12 (6.24) 

C3 = 0.3670959657 C4 = — tt^ - -. 

18 3 

Using equations ( |6.19 ) and ( 6.20 ), the new A± and o"^ are plotted as a function of 
u (figure [6.41 ). Complex solutions exist for uj < 0.595 where (as with the gaussian 
trial function) the discriminant in ( |6.19| ) is negative. The sech function also yields 
unphysical solutions for all of the "-" roots since a'^ is negative everywhere. Again, 
there is a part of the "+" solution (0.73 < a; < 1.41) that has positive al and can 
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a 



Figure 6.3: Plot of A± (top) and a\ (bottom) as a function of to for gaussian f{r/a). 
Ttie "+" roots are plotted in dashed green, while the "-" roots are plotted in solid red. 
The solution is not plotted for w < 0.437 as the discriminant in equation ( 3.19| ) is negative 
and only real A are of interest. The "+" roots are all unphysical since a^ is less than zero 
everywhere, which means that solutions within this region have an imaginary width. The 
solution for ct^ is also unphysical outside the region 0.62 < lj < 1.41 for the same reason. 
Fortunately, there are solutions that are both physical (real roots) and similar to those seen 
in the dynamic simulations (near the x's). 
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Parameter 


Dynamic 
Solution 


A+ = Ad 


a+ = an 


A 
a 

UJ 


0.62 
2.2 

1.4 


0.62 
2.5 
1.29 


0.72 
2.2 

1.27 



Table 6.2: Table of free parameter values {A, a, and uj) determined using a sech trial 
function (third and fourth columns) compared to the values estimated from the spherically 
symmetric dynamics (second column); the D subscript indicates the value was obtained from 
the dynamic solution. The third column contains results obtained using the trial function 
approach where A^ equals the amplitude of the n = 1 mode dynamic solution; the fourth 
column contains results where cr+ equals the width of the n = 1 mode dynamic solution. 
The data from the third and fourth columns occur around the x's in figure 6.4. 



therefore support physical solutions. Table S^ compares sample solutions (located 
near the x's in figure |6.4| ) to estimates obtained from the dynamic solution. 



6.2.3 Hyperbolic Secant (Squared) Trial Functions 

Lastly, in this subsection the field configuration is assumed to have a hyperbolic 
secant squared shape, 



^li 



sech I — 
a 



With the new choice of f{r/a), the integrals ( |6.16| ) are now reevaluated: 



(6.25) 



ci 



C3 



-vr 



45 
0.1053157513 



C2 



C4 



— vr 

18 3 

4 9 14 



(6.26) 



to 



105 45 

Using equations ( 6.19| ) and ( |6.20| ), the new A± and o"^ are plotted as a function of 
(figure |6.5D . Complex solutions exist for tu < 0.530 where the discriminant in ( |6.19| ) 
is negative. The sech function also yields unphysical solutions for all of the "-" 
roots since o"^ is negative everywhere. Again, there is a part of the "+" solution 
(0.69 < UJ < 1.41) that has positive ul and can therefore support physical solutions. 
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Plot of A± (top) and a\ (bottom) as a function of w for hyperbolic secant 
+" roots are plotted in dashed green, while the "-" roots are plotted in so lid 



Figure 6.4: 
/(r/a). The 

red. The solution is not plotted for u> < 0.595 as the discriminant in equation (6.19) is 
negative and only real A are of interest. The "+" roots are all unphysical since a'j^ is less 
than zero everywhere, which means that solutions within this region have an imaginary 
width. The solution for a^ is also unphysical outside the region 0.73 < a; < 1.41 for the 
same reason. Fortunately, there are solutions that are both physical (real roots) and similar 
to those seen in the dynamic simulations (near the x's). 
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Parameter 


Dynamic 
Solution 


A+ = Ad 


a+ = (TD 


A 
a 

UJ 


0.62 
3.4 

1.4 


0.62 
2.5 

1.28 


0.36 
3.4 

1.34 



Table 6.3: Table of free parameter values [A^ a, and uj) determined using a sech^ trial 
function (third and fourth columns) compared to the values estimated from the spherically 
symmetric dynamics (second column); the D subscript indicates the value was obtained from 
the dynamic solution. The third column contains results obtained using the trial function 
approach where 71+ equals the amplitude of the n ~ 1 mode dynamic solution; the fourth 
column contains results where a+ equals the width of the n = \ mode dynamic solution. 
The data from the third and fourth columns occur around the x's in figure p. 



Table O compares sample solutions (located near the x's in figure 6^) to estimates 
obtained from the dynamic solution. 



117 



15 

10 

^ 5 



< 




-5 

10 




CJ 




a 



Figure 6.5: Plot of A± (top) and a\ (bottom) as a function of w for hyperbolic secant 
squared f{r/a). The "+" roots are plotted in dashed green, while the "-" roots are plotted 



in solid red. The solution is not plotted for to < 0.530 as the discriminant in equation (6.19) 



is negative and only real A are of interest. The "+" roots are all unphysical since a^ is 
less than zero everywhere, which means that solutions within this region have an imaginary 
width. The solution for a^ is also unphysical outside the region 0.69 < a; < 1.41 for the 
same reason. Fortunately, there are solutions that are both physical (real roots) and similar 
to those seen in the dynamic simulations (near the x's). 
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Chapter 7 



Conclusions and Future Work 



Despite its simple appearance and long history, the real scalar field keeps revealing 
useful and fascinating physics to those who continue to explore it. Scalar field theory 
(both classical and quantum) is being used throughout the frontiers of physics; al- 
though this work stresses its application to cosmology, the phenomenology discussed 
is widely applicable to the general theory of phase transitions. 

Although oscillons had been discovered in 1977 0,[^, and analyzed in slightly 
more detail in 1995 [24|, the limitations in computation and numerical methods avail- 



able to both research teams prevented thorough exploration|^. This work presents a 
new technique to absorb outgoing radiation in a stable and non-reflecting manner 
that works well for massive flelds. This numerical method allows for the equations 
of motion to be evolved on a (relatively) small static lattice, making evolutions 
that previously required computational resources proportional to the square of the 

oscillon lifetime need only computational resources directly proportional to the oscil- 

^ Of course, limited computational resources was more of a factor for Bogolyugskii than Copeland 
et al., but the numerical methods that both were employing were inefficient. 
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Ion lifetime. This efficiency is achieved with the combination of solving the system 
of equations in MIB coordinates while including higher order dissipation. In MIB 
cooridinates outgoing (and ingoing) characteristic velocities go to zero near the 
outer edge of the computational domain so the outgoing waves are compressed to 
nearly the Nyquist limit on the lattice where they are then quenched by dissipation. 
The technique proves to be stable and free of contamination in both one and two 
dimensional problems. 



Following the more recent work of Copeland et al. |24], but with the aid of 
this new numerical method, this work reveals many new properties of oscillons (dis- 
cussed in chapter Q). In spherical symmetry and using the symmetric double- well 
potential, the method was integral to the discovery of a new resonant phenomenon. 
Given the long lifetimes of oscillons (particularly resonant oscillons) the MIB coor- 
dinate evolutions are as much as one hundred times more computationally efficient 
than those previously employed. The resonances display a time scaling law and 
reveal (on the threshold of successive shape mode modulations) a new type of non- 
radiative oscillon solution. This solution's existence and form is verified by solving 
a system of ordinary differential equations obtained from a periodic non-radiative 
ansatz. The method is also used to explore the threshold of expanding bubble for- 
mation for bubble collapse with an asymmetric double-well potential. A time scaling 
law is shown to exist with a scaling exponent, 7 ~ 2.1. 

Chapter |5| discusses the use of MIB coordinates in two dimensional axisym- 
metric bubble collisions. Evolutions are observed that result in expanding true- 
vacuum bubble formation, annihilation, coalescence, or soliton-like transmission. 
By looking at a two dimensional parameter space, (collision) boost velocity versus 
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bubble width, the threshold of expanding bubble formation is again explored. The 
threshold reveals a time scaling law with a scaling exponent (7 ~ 2.2) similar to the 
spherically symmetric collapse (7 ~ 2.1). This scaling law appears to be universal 
in that the scaling exponents observed are independent both of the values of the 
critical parameters {a* and v^) and which parameter is varied to explore the thresh- 
old (ie. T = jln\p — p*\ for p equal to a or vb)- The similarity between the scaling 
parameters of the spherically symmetric collapse and the axisymmetric collisions 
suggests that the dominant unstable mode (responsible for the critical phenomena 
observed) in the axisymmetric collisions might be spherically symmetric. In addi- 
tion to more critical phenomena, the parameter space surveys reveal that there are 
many collision-induced expanding bubbles. In other words, there are many bubble 
configurations that would otherwise have formed oscillons and eventually dispersed, 
but when colliding with other bubbles combine to form an expanding vacuum bub- 
ble. In a system where bubble collisions might be occurring, nucleation rates for 
expanding bubbles will be higher than those originally believed to arise solely from 
thermal (or quantum) fluctuations. 

Chapter g discusses oscillon solutions obtained using trial function methods. 
Judicious choices for trial functions allow the complete integration of the action; 
resulting in a system of algebraic equations instead of differential equations. Al- 
though the solutions obtained are not solutions to the true equations of motion, the 
predicted values for the amplitude, frequency, and width of the resonant oscillons 
are in reasonable agreement with observation. Furthermore, since the non-radiative 
resonant solution discussed in chapter was not known before this work, the mere 
existence of the non-radiative oscillon solution obtained in such a straightforward 
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manner was exciting. 

Although there is certainly more to be discovered within the nonlinear Klein- 
Gordon equation, it would be interesting to see if similar physics can be found in 
other models. It would be interesting to see how the stability of the oscillon solu- 
tions are affected by the coupling to gravity (GR). The extra gravitational attrac- 
tion might make oscillons more stable, however, the extra nonlinearities introduced 
might have entirely the opposite effect. It would also be interesting to see if similar 
phenomena are observable in the charged scalar field, both alone and coupled to elec- 
tromagnetism (the charged scalar field with a SDWP coupled to electromagnetism 
is a simple model of superconductivity). Lastly, the numerical methods used here 
might be applicable to the solution of other nonlinear partial differential equations. 
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Appendix A 



Will Cosmological Oscillons 
Form Black Holes? 



One of the most interesting applications of oscillons is to cosmological phase tran- 
sitions resulting from spontaneously broken gauge symmetries. But whenever one 
discusses the dynamics of massive objects in a cosmological context, two questions 
naturally arise: How large are these objects? and How dense are these objects?. If 
the object is cosmologically large, then gravitational effects might need to be in- 
cluded since the expansion of the universe is measurable on the length-scale of the 
object being studied. If the object is very dense, then there will be gravitational 
self-interaction that will likely change the dynamics of the system. This appendix 
answers these two questions and shows that oscillons do not need to be modeled 
with general relativistic effects included. 

The first thing to remember in this calculation is that (in contradiction to 
most relativist's conventions), this thesis uses typical high energy physics units where 
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h = c = 1, which imphes that energy can be treated as the fundamental dimension: 

[Energy] = [Mass] = [Length]"^ = [Time]"^ (A.l) 

Two useful conversion factors are IGeV = 1.7827 x lO"^^^ and IGeV'^ = 1.9733 x 

10~^^cm. 

This thesis almost exclusively uses dimensionless units, obtained by taking 

r = prriH, t = rniH, and (j) = V"? where mn is the mass in the model (the 

iriH 

Higgs mass) and p, t, and ^ are the dimensionful length, time, and scalar field, 

respectively]^. Converting back to dimensionful quantities, the typical oscillon radius 

and mass/energy are r ps 2.bm'^ and M ss 43——, respectively. Assuming the Higgs 

A 

mass is known in GeV, converting the oscillon radius radius to cm is performed by 

r = 2.5'mjj 

2.5 \ /l.97x 10-^^GeVcm\ 

^) [ 1 ) (^-2) 

'4.93 X IQ-'^'^GeV cm\ 



\ ruH J 

so that for an "electroweak" oscillon where the Higgs mass is m^ ~ lOOGe^, r ~ 
4.93 X lO^^^cm, or for a GUT oscillon with mn « lO^^GeV, r « 4.93 x 10"^^ cm. It 
is clear that these are not cosmological in their size, but what about their densities? 
Since these "particles" are so small, it might be inappropriate to even con- 
sider their classical gravitational Schwarzchild radius, but for an oscillon with M = 
43mH (taking A = 1), the electroweak oscillons have M ~ 8 x lO^^"'^^ while a GUT 

oscillon would have M ~ 8 x 10"^^. These oscillons then would have Schwarzchild 

-50 ^_ ( ^h\ _ ^^ _ in-48^_ „„j ^„ _ in-35. 



radii of rg = 1.13 x 10 cm , or r5 ~ 10 cm and rg ~ 10 cm, respec- 

\GeV , 

tively. 



^Also remember A is always dimensionless. 
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Oscillons 


Formula 




Attribute 


(Dependence on 


itih) 


Radius 


r = 4.9 X IQ-^'^cm 


fGeV\ 


Lifetime 


r = 3.3x lO^^^s 


'GeV\ 
K ruH ) 


Mass 


M = 7.7 X lO'^^g 


\Gev) 


Schwarzschild Radius 


rs = 1.1 X 10"^° cm 


\GeVj 



Table A.l: Table displaying how to determine typical oscillon attributes them as a function 
of the Higgs mass. The formulae were derived using typical values (in the dimensionless 
coordinates used throughout this thesis) of r = 2.5, T = 5000, and M ~ 43, for the oscillon 
radius, lifetime, and mass, respectively. 

So, although we do not know exactly how to describe all the forces of nature 
acting together at GUT energy scales, the point to be made here is that if the 
nonlinear Klein-Gordon action is minimally coupled to gravity, oscillons do not have 
a radius comparable to their Schwarzchild radius until the Higgs mass approaches 
ruH ^ W^^GeV. 
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